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ABSTRACT 


This  thesis  comparer  the  detection  perfonnance  of  the  1-1/2  D  instantaneous  power 
spectrum  (1-1/2  Djpj),  the  bispectrum,  the  instantaneous  higher-order  moment  slice 
(IHOMS)  method,  and  the  spectrogram  for  multi-component  stationary  signals, 
harmonically  related  stationary  signals,  and  multi-component  linear  FM  signals  corrupted 
by  additive  white  Gaussian  noise.  In  addition,  a  determination  of  the  relative  processing 
gain  between  the  1-1/2  Djpg  method  and  the  spectrogram  is  made  for  stationary  signals  in 
noise. 

The  results  of  this  thesis  show  that  1-1^  Djp^  has  a  processing  gain  advantage  over 
that  of  the  spectrogram  for  a  range  of  input  SNR  that  depends  upon  the  size  of  the  data 
window.  Under  some  conditions,  the  bispectrum  can  detect  both  harmonic  coupling  and 
phase  coupling  between  the  components  of  multi-component  signals.  IHOMS’  ability  to 
detect  linear  chirps  in  noise  is  limited  to  chirps  having  different  slew  rates,  and  the  method 
has  a  significantly  greater  computational  cost  than  both  the  spectrogram  and  1-1/2  Djpg. 
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I.  INTRODUCTION 


A.  OVERVIEW 

Spectral  estimation  is  an  important  data  analysis  tool  used  in  signal  processing  to 
determine  the  disTibution  of  power  as  a  function  of  frequency.  The  “classical”  Fourier 
transform  methods  and  a  variety  of  more  recently  developed  parametric  methods  are  now 
commonly  used  to  estimate  power  spectral  density.  These  methods  are  based  on  the 
autocorrelation  domain  and  therefore  rely  upon  information  contained  in  the  second  order 
moment  of  the  data  sequence.  Recent  improvements  in  computational  capability  have 
opened  up  a  new  and  active  research  area  concerned  with  the  extraction  of  additional 
information  contained  in  the  data  sequence’s  higher-than-second  order  statistics.  [Refs. 
1.2] 

Higher  order  statistics  involve  cumulants  rather  than  moments.  Cumulants  and 
moments  are  similar,  and  each  can  be  expressed  in  terms  of  the  other.  In  fact,  under 
certain  conditions  they  are  identical.  Higher  order  spectra  (HOS),  also  known  as 
polyspectra,  are  formed  by  taking  the  Fourier  transform  of  the  cumulant  sequence.  For 
example,  the  two-dimensional  Fourier  transform  of  the  third  OTder  cumulant  yields  the 
polyspectrum  known  as  the  bispectrum.  Likewise,  the  three-dimensional  Fourier 
transform  of  the  fourth  order  cumulant  is  called  the  trispectrum.  If  the  first  moment  is 
zero,  the  second  order  cumulant  is  equal  to  the  second  order  moment,  and  the  one¬ 
dimensional  Fourier  transform  of  either  sequence  produces  the  familiar  power  spectrum. 
[Refs.  1,3] 

Higher  order  statistics  and  their  spectra  possess  some  exploitable  properties  that  their 
second  order  counterparts  do  not.  One  potential  advantage  involves  a  greater  resistance 
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to  the  effects  of  certain  types  of  additive  noise.  For  white  or  colored,  zero  mean  Gaussian 
processes,  the  cumulant  of  any  order  greater  than  two  is  zero  and  the  corresponding 
polyspectrum  is  zero.  Consequently,  polyspectra  are  expected  to  be  more  resistant  than 
the  power  spectrum  is  to  the  effects  of  additive  Gaussian  noise.  One  limitation  of  the 
power  spectrum  is  that  phase  cannot  be  accurately  reconstructed  unless  the  underlying 
signal  or  system  is  minimum  phase.  In  contrast,  higher  order  spectra  preserve  true  phase 
information.  Finally,  HOS  are  useful  for  detection  and  classification  of  non- 
linearities.  Depending  upon  the  application,  the  extra  computational  cost  of  HOS  may  be 
justified.  [Refs.  3,4] 

The  bispectrum,  the  1-1/2  D  instantaneous  power  spectrum  [Ref  5],  and  the 
mstantaneous  higher  order  moment  slice  (IHOMS)  method  [Ref.  6]  are  studied  in  this 
thesis.  Performance  comparisons  are  made  with  respect  to  each  method  as  well  as  to  an 
accepted  second  order  method. 

B.  THESIS  OUTLINE 

The  following  describes  the  organization  of  the  remainder  of  this  thesis.  Chapter  n 
first  defines  cumulants  and  shows  their  relation  to  moments.  The  essentials  are  then 
presented  for  each  spectral  estimation  method  studied.  Chapter  m  makes  a  processing 
gain  comparison  between  the  1-1/2  D  instantaneous  power  spectrum  and  the  spectrogram. 
Chapter  FV  displays  simulation  results  that  show  how  the  different  methods  perform  on  a 
variety  of  signals.  Chapter  V  presents  conclusions  and  suggestions  for  future 
work.  Appendices  contain  the  detailed  steps  of  Chapter  HI  calculations  and  the  Matlab 
programs  used  to  execute  the  computer  simulations. 


II.  ELEMENTS  OF  HIGHER  ORDER  SPECTRAL  ANLYSIS 


This  chapter  describes  the  higher  order  based  techniques  compared  later  in  this 
thesis.  S ufficient  information  is  presented  in  order  to  develop  an  adequate  understanding 
of  each  method.  We  first  define  cumulants  and  show  how  they  are  used  instead  of 
moments  to  form  polyspectra. 

A.  CUMULANTS  AND  MOMENTS 

In  general,  the  characteristic  function.^  0a>).  is  defined  as  the  conjugated  Fourier 
transform  of  a  random  variable’s  probability  density  function,  f(jc)  [Ref.  7]: 

oo 

^(M)  =  jfix)exp(ja>x)dx.  (2.1) 

The  right-hand  side  of  (2.1)  is  simply  the  expectation  of  exp  (JaX) .  The  n^  moment  of 

the  random  variable  X  can  be  found  by  evaluating  the  n**^  derivative  of  the  characteristic 
function  with  ©  set  equal  to  zero  [Ref.  7]: 

=  (-j)"^4>0©)  .  (2.2) 

0,-0 

In  comparison,  the  n*  cumulant  is  defined  as  the  nth  derivative  of  the  natural  logarithm  of 
the  characteristic  function  evaluated  at  ©  =  0  [Ref.  4]: 

C{X"}  =  (-;)”^n[<t>(®)]  (2.3) 

<'“>  „.o 

Although  (2.3)  clearly  shows  the  basic  difference  between  a  moment  and  a  cumulant, 
it  does  not  readily  show  that  a  cumulant  can  be  considered  as  a  moment  with  its  lower  order 
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statistical  dependence  removed  [Ref.  8].  This  is  perhaps  best  seen  by  computing 
cumulants  in  an  alternate  manner  that  involves  a  partitioning  of  equal  and  lower  order 
moments. 

Given  a  set  of  random  variables  X  =  Xj, ....  ,  Ix  =  {1.2 . k}  is 

defined  to  be  the  set  of  component  indices  contained  in  X.  Denoting  a  subset  of  7;^  by  /, 
using  nixU)  to  represent  the  moment  of  those  components  of  X  comprising  Xj,  and 
denoting  the  cumulant  of  Xj  as  CxU) .  the  moment-to-cumulant  (M-C)  equation  can  be 
written  as  [Ref.  3]; 

=  L  (2.4) 

where  Ip  are  non-intersecting  non-empty  subsets  of  7  that  form  the  partitions,  and  q  is  var¬ 
ied  from  one  to  the  number  of  elements  in  7.  The  summation  in  (2.4)  is  over  all  unique 
partitions  and  for  all  q.  The  first  order  cumulant  is  equal  to  the  first  order  moment  since 
only  one  unique  partition  exists  for  one  random  variable.  Table  2.1  shows  how  (2.4)  is 
applied  to  compute  the  moment  partitions  that  form  the  second  order  cumulant.  Table 
2.2  displays  the  details  of  the  third  order  cumulant  calculation.  Summing  the  last  column 
of  Table  2.1  yields  the  second  order  cumulant.  Summing  the  last  column  of  Table  2.2 
produces  the  third  order  cumulant.  The  first  three  orders  of  cumulants,  expressed  in 
terms  of  the  component  indices  of  a  set  of  random  variables,  are  then: 
c(l)  =  mil), 

c(l,2)  =  m(l,2)-m(l)m(2) , 

c(l,2,3)  =  m(l,2,3)  +2m(l)/n(2)m(3)  -m(l,2)m(3) 
-m(l,3)m(2)  -m(2. 3)m(l). 

Calculation  of  the  fourth  order  ciunulant  is  similarly  detailed  in  [Ref.  3].  An  expression 
that  calculates  moments  from  partitioned  cumulants  is  also  given  in  [Ref.  3]. 
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TABLE  2.1:  COMPUTATION  OF  SECOND  ORDER  CUMULANT 


h 

h 

q 

M-C 

equation 

(2.4) 

1.2 

1 

m(1.2) 

1 

2 

2 

-m(l)m(2) 

TABLE  2.2;  COMPUTATION  OF  THIRD  ORDER  CUMULANT 


h 

h 

q 

M-C  equation  (2.4) 

1.2.3 

1 

m(  1,2,3) 

1.2 

3 

2 

-m(l,2)m(3) 

1.3 

2 

2 

-m(l,3)m(2) 

2.3 

1 

2 

-m(2,3)m(l) 

1 

2 

3 

3 

2m(l)m(2)m(3) 

Inspection  of  the  moment-cumulant  expressions  show  that  dependencies  among 
random  variables  are  removed  when  cumulants  are  calculated.  In  fact,  both  the  second 
and  third  order  cumulants  are  zero  if  all  the  random  variables  are  independent  of  each  other. 
[Ref.  8] 

Another  observation  regarding  the  above  moment-cumulant  expressions  applies  in  the 
commonly  encountered  situation  where  the  means  of  the  random  variables  are  equal  to 
zero.  When  this  is  true,  the  first  order  cumulant  is  zero  since  it  equals  the  mean,  the 
second  order  cumulant  simply  equals  the  variance,  and  the  third  order  cumulant  equals  the 
third  order  moment.  [Ref.  3] 
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There  is  more  than  one  acceptable  choice  when  it  comes  to  choosing  which  term(s)  to 
conjugate  in  the  expectation  expression  for  moments  of  order  greater  than  two.  Different 
choices  have  different  consequences  when  the  signal  being  processed  is  complex  [Refs.  9. 
10].  One  consequence  is  considered  later  when  the  symmetry  region  of  the  bispectral 
plane  is  discussed.  For  now,  the  conjugated  terms  are  chosen  arbitrarily.  The  following 
expressions  show  the  details  of  the  moment  calculations  in  the  last  set  of  equations  for  the 
zero  mean  case  [Ref.  1 1]: 

=  £{jc[n]}  =  0,  (2.5) 

=  E{x*[n]x[n  +  l^]},  (2.6) 

=  E{x*[n]x[n  +  li]x[n  +  l2]}.  (2.7) 

The  expression  for  the  fourth  order  cumulant  as  calculated  by  (2.4)  has  15  moment  terms. 
If  the  mean  is  zero  all  but  four  terms  go  to  zero  when  the  signal  is  real  [Refs.  3,11]: 

[^1.  h  =  E{x[n]x[n  +  li]x[n  +  l2]x[n’i-  /j] } 

-E{x[n]xln  +  ll]}E{x[n  +  l2]x[n  +  l^] } 

-£  {X  [n]  jc  [n  +  /j] }  £  {X  [n  +  /j]  X  [n  +  /j] }  (2.8) 

-£  {x  [n]  x[n  + 1^]}  E  {x[n  +  l^]x[n  + 12] } . 

Provided  that  two  of  the  four  random  variables  are  conjugated  in  the  expectation  terms 
of  (2.8),  the  fourth  order  cumulant  of  a  zero  mean  complex  process  has  just  three 
expectation  terms.  One  of  the  last  three  terms  in  (2.8)  will  be  zero.  The  zero  term 
depends  upon  which  variables  are  conjugated.  [Ref.  11] 

Cumulants  have  three  properties  that  make  them  more  desirable  than  moments  when 
it  comes  to  higher  order  statistics: 

1.  Each  cumulant  is  independent  of  all  lower  order  cumulants.  Consequently,  all 
cumulants  of  order  greater  than  two  are  equal  to  zero  for  a  Gaussian  process,  as 
a  Gaussian  process  is  completely  characterized  by  its  first  and  second 
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moments.  Higher  order  moments,  on  the  other  hand,  can  contain  information 
about  lower  order  moments.  In  fact,  only  higher  order  moments  of  odd  order 
are  zero  for  a  Gaussian  process.  The  non-zero  even  order  higher  moments  do 
not  contain  any  information  not  already  contained  in  the  first  two  moments. 
[Ref.  8] 

2.  If  a  set  of  n  random  variables  can  be  divided  into  more  than  one  statistically 
independent  subset,  then  the  n*^  order  cumulant,  unlike  the  n^  order  moment,  is 
equal  to  zero.  [Ref.  4] 

3.  Unlike  moments,  the  cumulant  of  the  sum  of  two  independent  stationary  random 
processes  is  equal  to  the  sum  of  the  cumulants  of  each  process.  [Ref.  4] 


B.  THE  BISPECTRUM 

1.  Deflnition  and  Region  of  Support 

The  bispectrum  is  defined  as  the  two-dimensional  Fourier  transform  of  the  third 
order  cumulant  [Ref  4]: 

oo  «0 

B(©j,©2)  =7^  L  +  •  (2-9) 

Analogous  to  the  power  spectrum,  the  bispectrum  can  also  be  defined  with  frequency 
domain  quantities  Given  N  samples  of  stationary  signal  r  ( n) ,  its  Fourier  transform  is 

N-l 

X(©)  =  Y^x{n)expi-join) .  (2.10) 

«  -  0 

Assuming  that  the  third  order  cumulant  in  (2.9)  is  computed  with  the  conjugation  scheme 
shown  in  (2.7),  the  equivalent  frequency  domain  expression  for  the  bispectrum  is  obtained 
through  an  extension  of  the  periodogram  [Refs.  9,12]: 
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(2.11) 


fi(o>j,C02)  =  ^X(0)j)X(C02)A*  ((0j+C02); 

where  |cOj|  ^  n,  |(02|  ^  n,  and  jtOj  +  ©21  ^  n.  Equations  (2.9)  and  (2.11)  represent  two 
different  non-parametric  methods  that  can  be  used  to  calculate  the  bispectrum.  The 
{q)proach  taken  by  (2.9)  is  called  the  indirect  method,  while  the  approach  used  by  (2.11)  is 
known  as  the  direct  method.  [Ref.  1] 

In  general,  the  bispectrum’s  region  of  support  is  a  hexagon  centered  at  the  origin 
of  the  (©j,  ©2)  bi-frequency  plane.  Evaluation  of  the  mean  product  of  three  Fourier 

amplitude  terms  shows  that  the  bispectral  plane  exhibits  certain  symmetries.  For  a 
stationary,  real,  continuous  time  signal  the  expected  value  of  the  product  of  three  Fourier 
components  is  [Ref.  13]: 

£{.Y^(©j)X^(©2)X^(©3)}  =  fi^(©j,©2)6(©j  +  ©2  +  ©3);  (2.12) 

where  the  subscript  c  denotes  continuous  time  quantities. 

Three  properties  of  (2.12)  determine  bispectral  symmetry.  First,  since  the 
frequency  indexes  must  sum  to  zero,  ©3  =  -  ©j  -  ©2 .  Second,  the  frequency  indexes  in 
the  expectation  operation  can  be  interchanged.  Third,  for  the  bispectrum  of  a  real  signal, 
conjugate  symmetry  results  since  B(-©j, -©2)  =£*(©^,©2).  Using  the  first 
property  to  express  ©3  in  terms  of  ©^  and  ©2,  the  symmetry  lines  for  a  continuous  time 
signal  are  [Ref.  13]; 
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II 

<• 

2©j  =  -©2 

(from  ©j 

= 

2©2  =  -©1 

(from  ©2 

II 

8 

3 

1 

II 

3 

0 

II 

3 

(from  ©j 

II 

1 

8 

8 


cOj  =  0 


(from  cOj  = 


2n 


If  the  signal  is  band-limited  to  to^  and  sampled  at  a  rate  of  (o^  =  2co^  =  — ,  the 
bispectrum  term  in  (2.12)  becomes  fi(©j  +  ^^,©2+ applying  the 


approximation; 


X((o)  s  £  x^(co+^) 

nm-eo 

2nn 


(2.13) 


The  frequency  indexes  in  (2.12)  now  must  sum  to  instead  of  zero.  As  a  result,  the 

following  set  of  symmetry  lines  are  formed  in  the  bispectrum  of  a  sampled  signal  [Ref. 
13]: 


OOj  =  C02* 


2nn 


2co,  +C0,  = 

1  2  f 


2nn 


+  2©2  =  ^ 


©1  =  -©2, 
2nn 

©2  =  — . 
2nn 


©,  =  - . 

t  T 

Figure  2.1  shows  the  region  of  support  and  its  symmetry  lines  for  the  first  of  an 
infinite  number  of  bispectral  hexagons  given  a  sampled  real  signal  [Refs.  12,13].  The 
sampling  period  is  assumed  to  be  equal  to  one.  A  complex  signal  has  the  same  hexagonal 
support  region  as  that  of  a  real  signal  but  it  is  symmetrical  about  only  one  of  three  dotted 
symmetry  lines  shown  in  Figure  2. 1 .  Depending  upon  the  term  conjugated  in  (2.7)  used  to 
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compute  the  third  order  cumulant  for  the  indirect  method,  or  which  term  in  (2.11)  is 
conjugated  to  compute  the  bispectrum  by  the  direct  method,  the  bispectrum  of  the  complex 
signal  exhibits  symmetry  in  only  one  quadrant  about  either  the  =  0)2  line,  the 

©2  =  -2©j  line,  or  the  =  -2©2  line.  The  conjugation  scheme  used  exclusively  in 
this  study  is  shown  in  (2.7)  and  (2.1 1).  Under  this  scheme,  the  bispectrum  of  a  complex 
signal  displays  two-fold  symmetry  in  the  first  quadrant  about  the  ©^  =  ©2  line.  There 

are  two  other  possible  ways  to  conjugate  just  one  term,  and  three  possible  ways  to  conjugate 
any  two  terms  in  the  applicable  expressions.  The  bispectrum  symmetry  for  these  other 
schemes  are  considered  fully  in  [Ref.  9}. 
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Because  of  the  symmetry  described  above,  only  a  small  portion  of  the  bispectrum 
needs  to  be  computed.  Once  the  bispectrum  is  known  for  this  small  region  known  as  the 
non-redundant  region,  the  rest  of  the  bispectrum  can  be  found  using  symmetry.  The  non- 
redundant  region,  again  for  the  sampled  signal  case,  is  shown  in  Figure  2.2.  [Ref.  12.13] 


Figure  2.2:  The  non-redundant  region  of  the  bispectrum. 

2.  Computation  by  the  Indirect  Method 

The  indirect  method  requires  a  two-dimensional  Fourier  transform  of  an 
estimated  third  order  cumulant  sequence.  A  »wo-dimensional  window  can  be  used 
provided  that  it  satisfies  certain  requirements  based  on  cumulant  symmetries.  The 
following  procedure  outlines  the  indirect  method  given  the  data  set 
{X(0),X(l),...X(A^-l)}:[Ref.  1] 

1 .  Form  K  segments  of  M  samples  such  that  N=KM. 

2.  Remove  the  mean  from  each  segment. 

3.  Defining  i=1.2....K  to  indicate  the  segment  number,  estimate  the  third  order 
moment  sequence  for  each  segment,  {j:  (it) ,  =  0, 1, . . .  Af  -  1 } ; 
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(2.14) 


£  [jc^‘^(/)]  il  +  m)x^‘'  (l  +  n); 


.(0 


(0 


M 


l  •*  s. 


where  sl=max(0,-m,-n)  and  s2=min(M-l,M-l-m,  M-l-n).  This  sequence  is 
also  called  the  tri-correlation  sequence. 

4.  Average  all  segment  moment  estimates; 

tHm,n)  =  p  V  (m.  n)  .  (2.15) 

5.  Multiply  the  averaged  moment  estimate  by  a  two-dimensional  window,  W(m,n), 
and  take  the  two-dimensional  Fourier  transform  to  obtain  the  bispectrum  esti¬ 
mate; 


L  L 

~  ^  ^(/n, n)W''(m, /i)exp{-7(Q)jm  +  a)2”) } ;  (2.16) 

mm~L  nm-L 

where  L<M-\. 

The  window,  W{m,n),  used  in  the  tri-correlation  sequence  defined  in  (2.16)  must 
satisfy  the  following  conditions  [Ref.  1,11]; 

1 .  The  window  conforms  to  third  order  cumulant  symmetry.  Specifically, 

W{m,Ti)  =  W{n,m)  =  W{-m,n-m)  =  W{m-n,-n). 

2.  The  window  is  zero  outside  the  region  of  support  for  the  tri-correlation 
sequence. 

3.  The  window  is  normalized  to  one  at  m=n=0. 

4.  The  Fourier  transform  of  the  window  is  non-negative. 

Windows  used  in  this  study  belong  to  a  separable  class  of  window  functions  defined  by 

W{m,n)  =  dim)d(n)din-m) ,  (2.17) 

where  an  appropriate  one-dimensional  window,  d(m),  was  chosen  subject  to  the  following 
constraints; 
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d(m)=^d(-m) 
d(m)=0  for  m>L, 
d(Ohl, 

D{(i))  ^0  for  all  co.  [Ref.  11] 

The  unit  hexagonal  window  was  constructed  by  applying  (2.17)  to  a  one¬ 
dimensional  rectangular  window  [Ref.  14], 

d (m)  =s  1  for  \m\  ^L, 
d(m)  =0  otherwise. 

The  one-dimensional  Parzen  window,  given  by 

2  3 

d{m)  =  l-6(-^)  +6(^)  forimi^2* 

3 

d{m)  =  2(1-^)  for^^\m\^L, 

d  (m)  =  0  for  \m\  >L, 

was  similarly  transformed  into  a  two-dimensional  Parzen  window  [Ref.  1].  The  last 
separable  window  considered  is  known  as  the  minimum  bias  supremum  window,  or 
simply  the  optimum  window.  The  corresponding  one-dimensional  window  is  defined  by 
[Ref.  1]: 


d(m)  =  i 

d(m)  =  0 


•  /"'"a 

Sin(-j-) 


for  \m\  >L. 


for  \m\  ^L, 


A  trade-off  situation  exists  between  the  Parzen  window  and  the  optimum  window 
with  regard  to  the  windows’  statistical  properties.  The  bias  of  the  optimum  window  is 
about  18%  smaller  than  that  of  the  Parzen  window.  However,  the  optimum  window  has 
a  26%  greater  variance.  [Ref.  15] 
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3.  Computation  by  the  Direct  Method 

The  direct  method  generates  the  bispectrum  hrom  (2.11).  This  requires  a  one- 
dimensional  Fourier  transform  of  the  signal.  Frequency  domain  smoothing  can  be  used 
to  reduce  the  variance  of  the  final  bispectrum  estimate.  The  direct  method  is  outlined  for 
the  data  set  {X(0),X(1).  1))  [Ref.  1,12]: 


1.  The  data  sequence  is  segmented  into  K  segments.  The  mean  is  removed  from 
each  segment  before  zero  padding  the  segment  to  a  convenient  fast  Fourier 
transform  length,  L. 

2.  An  L-point  Fourier  transform,  (X) ,  is  computed  for  each  segment.  The 
superscript  designates  the  segment  number. 


3.  The  bispectrum  estimate  of  the  segment,  denoted  as  b,  is  then  computed  by: 

5.(Xi,X2)  = 

4.  Two  different  approaches  can  be  taken  now  to  reduce  the  variance  of  the  final 
estimate.  The  first  approach,  denoted  Bi,  involves  smoothing  each  segment 
separately  and  then  averaging  the  smoothed  segment  estimates.  The 
smoothing  is  accomplished  uniformly  over  each  bispectral  location’s 
neighboring  frequency  points  [Ref.  12]: 

R/2  R/2 

5i(©p®2)  ^  L  L  (^1  +  ^ ^2 (2.18) 

^  /  -I  «  rm-R/2  s^-R/2 


where  ©j  =  (-^)Xj,©2  =  )X2,  and/?^L. 

In  the  second  {q>proach,  denoted  B2,  all  the  unsmoothed  segment  estimates  are 
first  coherently  averaged; 


82  (©p  ©2) 


1  ~ 

^£^(©^,©2). 


(2.19) 
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Then  a  frequency  domain  smoothing  window  is  applied  through  a  two- 
dimensional  convolution.  [Ref.  14] 

Hi-Spec  software  [Ref.  14]  is  used  in  this  study  to  compute  the  bispectrum  by  the 
direct  method.  This  program  uses  the  second  approach  described  above,  based  upon 
(2.19)  to  reduce  the  variance  of  the  fmal  bispectrum  estimate.  There  are  a  few  cations 
available  in  this  software  package  to  accomplish  the  frequency  domain  smoothing.  The 
option  chosen  most  frequently  is  the  Rao-Gabr  window  defined  by: 


W  (m,  n) 


m^  +  n^  +  mn 


for  (m,n)  £  G; 


where  M  =  2  •  Fourier  transform  length.  The  set  G  is  the  collection  of  (m,  n) 

points  satisfying 


1  7 

m  +n  +mn^ 


winlen 


The  parameter  winlen  is  the  desired  window  length.  Spatial  smoothing  is  applied  over 
the  (winlen)  ^  adjacent  frequency  points  around  each  bispectrum  point.  [Refs.  12,14] 


4.  Comparison  of  Methods 

Both  the  direct  method  and  the  indir^t  method  lead  to  asymptotically  unbiased 
estimators,  i.e.,  E  {B  (o^,  ©2) }  ^  fi  (©j,  ©2)  .  The  indirect  method  has  an  asymptotic 
variance  given  by: 

var  {Re  B  (©j,  ©2) }  =  var  {Im  B(©j,  ©2) } 

-  — ^F(©j)F(©2)F(©j+©2); 

M  K 

where  is  the  energy  of  the  tri-correlation  domain  window,  F(©)  is  the  true  power 
spectrum,  Af  =  2L  +  1  is  the  Fourier  transform  size  (i.e.,  an  MxM  FFT),  and  AT  is  the 
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number  of  segments  averaged  [Ref.  1].  The  asymptotic  variance  of  the  direct  method 
when  using  the  smoothing  approach  described  by  (2.18)  is; 

var  {Re  (coj,  CD2) }  =  var  (Im  (tOj,  002) } 

=  -^/*(a)i)/"(cD2)^(“i  +  o)2); 

2KL 

where  is  the  length  of  the  data  sequence,  K  is  the  number  of  segments,  and  L  is  the  size 
of  each  segment.  [Ref.  12] 

The  direct  method  and  the  indirect  method  are  identical  when  neither  uses  a 
window,  and  the  indirect  method’s  tri-correlation  estimate  is  computed  for  a  number  of  lags 
equal  to  the  full  segment  length  so  that  L  =  M-  \  in  (2.16). 

C.  THE  1-1/2  D  INSTANTANEOUS  POWER  SPECTRUM 

The  1-1/2  D  instantaneous  power  spectrum  (1-1/2  Djps)  is  a  combination  of  the 
standard  1-1/2  D  spectrum  (1-1/2  Dsyj)  and  the  Instantaneous  Power  Spectrum  (IPS)  [Refs. 
3,5,16].  It  is  shown  in  [Ref.  5]  that  1-1/2  Djpj  performs  better  than  the  conventional 
spectrogram  in  some  respects.  When  used  to  process  dynamic  signals,  1-1/2  Djpj  is 
observed  to  have  a  quicker  spectral  rise  time  and  a  quicker  fall-off  time  than  the 
spectrogram.  The  1-1/2  Djps  method  also  appears  to  be  good  at  detecting  low  SNR 
stationary  signals.  [Ref.  5] 

The  standard  1-1/2  D  spectrum  is  a  degenerate  form  of  the  bispectnim.  When  the 
first  lag  in  (2.7)  to  set  to  zero,  the  third  order  cumulant  expression  becomes 

C^^^(0./2)  =  £{^(n)X(n)X(n  +  /2)}  (2.20) 

=  E{\X{n)\h{n  +  l2)]. 

A  one-dimensional  Fourier  transform  of  (2.20)  produces  the  1-1/2  Dguj  spectrum.  [Ref.  3] 
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The  Instantaneous  Power  Spectrum  is  defined  as  the  average  of  the  derivatives  of 


Page’s  miming  spectrum. 


ou.n  =l 


I 

J  5  (x)  exp  (-j2nfx)  dx\ 


(2.21) 


and  Levin’s  backward  running  spectrum. 


=4 


p  (t)  exp  i-j2nfx)  dx 


(2.22) 


A  discrete  version  of  IPS  is  obtained  by  taking  the  Fourier  transform  of  windowed 
correlation  estimates  formed  using  simple  lag  products  vice  a  sum  of  lagged  products: 

N-\ 

/P5(/i.(d)  =  5  £  {^(n);^  (n-/t) +A*(/i)X(/i  +  /t)}  w(^■)e:tp(-y•o)^);  (2.23) 

where  w{k)  is  the  window  function,  and  N  is  the  length  of  the  sampled  data  sequence 
{X  (0) ,  X ( 1) , X (N  -  1) }  [Refs.  5,16]. 

The  1-1/2  Dips  algorithm  is  created  by  substituting  the  estimate  of  (2.20)  into  (2.23)  so 
that 

N-\ 

1-1/2  DjpJn,  ©)  =.£  {\X{n)\^2t  (n-k)  ■^\X{n)\h(n-^k)}w{k)exp{-jis,k) 

(2.24) 

with  variables  defined  as  in  the  IPS  equation.  [Ref.  5] 


D.  INSTANTANEOUS  HIGHER-ORDER  MOMENT  SLICE 

The  time-frequency  representation  introduced  in  [Ref.  6]  is  based  upon  a  one¬ 
dimensional  Fourier  transform  of  a  higher-order  moment  slice.  While  satisfying  certain 
constraints,  a  desired  number  of  different  moment  slice  sets  are  generated  in  order  to  form 
an  equal  number  of  different  but  related  time-frequency  surfaces.  As  a  consequence  of 
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the  moment  slice  construction  constraints,  signal  dependent  auto-terms  are  enhanced  and 
coss-terms  are  de-emphasized  when  the  tune-frequency  surfaces  are  coherently 
summed.  To  further  control  the  contrast  between  auto-terms  and  cross-terms,  each 
moment  slice  is  multiplied  by  a  kernel  function  before  being  Fourier  transformed.  The 
kernel  function  controls  the  degree  of  auto-term  and  cross-term  concentration  on  the  time- 
frequency  stuface.  [Ref.  6] 

The  instantaneous  higher-order  moment  slice  (IHOMS)  method  is  used  to  detect 
rnulticomponent  linear  chirp  signals  having  different  instantaneous  frequency  slopes,  or 
slew  rates.  One  characteristic  of  this  time-frequency  representation  is  that  auto-term  lines 
pass  through  the  origin  of  the  time-frequency  plane.  The  number  of  auto-terms  so 
identified  indicate  the  munber  of  chirps  in  the  signal.  In  addition,  the  angle  that  the 
auto-term  line  makes  with  the  time  axis  is  related  to  the  instantaneous  frequency  of  the 
chirp.  [Ref.  6] 

The  IHOMS  method  is  outlined  for  a  received  signal,  x  (r) ,  given  by 

M 

x{t)  =  £eJrp[/(M,.^-^h/)],  (2.25) 

<  - 1 

where  M  is  the  number  of  chirps  in  the  received  signal.  The  instantaneous  frequency  of 
the  i***  chirp  is  given  by: 

fiM  = 

A  general  expression  for  a  sampled  version  of  (2.25)  is: 

M  r  2  1 

x{n)  =  J^exp{j2n  (^)m,+ (^)  h,  };  (2.27) 
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where  the  sampling  period  is  assumed  to  be  unity,  and  N  is  the  number  of  discrete  samples 
in  the  data  record.  [Ref.  6] 

The  following  procedure  is  used  to  form  the  IHOMS  time-frequency  surface  [Ref.  6]; 


1.  The  desired  wder  of  the  moment  slice,  r,  is  chosen  first.  The  order  should  be 
greater  than  three  so  that  a  variety  of  related  time-frequency  surfaces  can  be 
sumxned  in  order  to  form  the  composite  surface.  In  general,  auto-terms  stand 
out  better  when  more  surfaces  are  summed.  However,  if  the  order  is  chosen 
too  large,  cross-term  effects  increase  and  may  cause  spurious  peaks  in  the  com¬ 
posite  time-firequency  surface. 


2.  The  parameter  vector,  g  = 


v-i 


,  is  generated  somewhat  arbi¬ 


trarily  but  subject  to  the  following  constraints: 


f-i 


I 

«- 1 


and 


(2.28) 


£  =  1. 


(2.29) 


3.  The  moment  slice  is  then  calculated  for  a  given  instant  of  time  (0,  a  given 
parameter  vector,  and  a  desired  number  of  lags  (x): 


m^('c;r,a^*^)  =  /  (T)x(/-hajT)  ...j:(t-ha^_jx) .  (2.30) 

4.  The  parameter  g  in  the  kernel  function,  exp  (-gx^) ,  is  chosen  very  small  so  that 
auto-terms  concentrate  along  straight  lines.  However,  if  g  is  made  too  small, 
cross-terms  show  up  everywhere  on  the  time-frequency  plane  and  tend  to  add  up 
like  auto-terms  when  the  surfaces  are  later  summed  to  form  the  composite 
surface. 
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5.  The  Fourier  transform  of  the  moment  slice  is  computed  and  multiplied  by  the 
kernel  function. 

oo 

Z  U,f,a )  =  J  exp  mj? (x;/,  )  exp  i-j2nfx)  dx.  (2.31) 

6.  Steps  three  through  five  are  repeated  for  a  different  time  instant,  the  same 
parameter  vector,  the  same  g,  and  the  same  number  of  lags.  Stacking  the 
results  of  (2.31)  for  successive  time  instants  forms  a  complex  valued  time-fre¬ 
quency  surface  based  upon  a  particular  parameter  vector. 

7.  A  new  parameter  vector  is  generated  and  then  steps  three  through  six  are 
repeated  to  create  another  time-frequency  surface.  This  is  done  again  and 
again  until  the  desired  number  of  complex  valued  surfaces  are  computed.  The 
composite  surface  is  obtained  by  taking  the  magnitude  of  the  coherently 
summed  parameter  vector  dependent  time-frequency  surfaces, 

K 

G(t,f)  =  £Z(t./;a^*))  .  (2.32) 

kml 

where  K  is  the  number  of  different  parameter  vectors  used  to  generate  an  equal 
number  of  related  time-frequency  surfaces. 

The  maxima  of  the  function, 

D(0)  =  ^  G(vcos0,  vsin0) ,  (2.33) 

vtlf, 

where  0  is  varied  between  zero  and  nil  radians,  and  the  radial  distance  v  is  varied 
incrementally  between  Ri  and  /?2.  correspond  to  chirps.  The  value  of  h  in  the  slope  term 
of  a  given  chirp’s  instantaneous  frequency  expression,  (2.26),  can  be  found  from  the  0 
producing  the  associated  maximum  value  in  (2.33)  using  the  relation: 
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where  the  subscript  j  denotes  a  particular  chirp.  ^  is  the  number  of  lags  used  to  compute 
the  moment  slice  in  (2.30),  and  N  is  defined  as  in  (2.27). 


III.  PROCESSING  GAIN  COMPARISONS 

In  this  chapter  we  compare  the  performance  of  1-1/2  Djps  to  that  of  the  periodogram 
and  its  time-frequency  representation,  the  spectrogram.  The  next  chapter  will  examine 
bispectral  as  well  as  IHOMS  methods,  and  compare  them  to  1- 1/2  Djp^  and  the  spectrogram 
where  qjpropriate. 

Matlab  programs  used  to  conduct  the  computer  simulations  for  this  cht^ter  are 
described  frrst.  Theoretical  calculations  and  simulation  results  are  then  presented  to  show 
how  each  method  treats  noise  and  signal  components  separately.  Finally,  relative  signal- 
to-noise  ratio  processing  gains  are  estimated. 


A.  SIMULATION  PROGRAMS 

The  1-1/2  Djpj  algorithm  was  originally  implemented  in  [Ref.  5]  with  the  extrinsic 
Matlab  function  ONE_HALF.  Modifications  were  subsequently  made  to  the  program  to 
remove  the  mean  from  each  data  segment  and  to  delete  some  unnecessary  smoothing 
steps .  The  revised  ONE_H ALF  function  is  contained  in  Appendix  C .  The  user  specifies 
an  input  data  sequence,  the  size  of  the  data  window,  a  step  increment,  and  a  window 
function.  ONE_HALF  returns  the  time-frequency  representation  of  the  input 
sequence.  The  ONE_HALF  algorithm  is  based  upon  (2.24)  but  is  constructed  in  such  a 
way  as  to  avoid  undesirable  product  term  interferences  that  may  form  spurious  peaks  in  the 
time-frequency  representation.  The  following  sequence  is  frrst  generated: 


prodl  = 


element{0)  element(l)  ...element -  1) 


where. 


(3.1) 
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element{k)  =  |x(/i)l^{jt*  (n-^) +jr(/i  +  it) } ,  (3.2) 

winlen  is  the  specified  data  window  length,  and  x(n)  is  the  input  data  segment.  A  second 
sequence,  prodl,  is  formed  by  deleting  the  first  element  in  prod\  and  then  reversing  the 
order  of  the  conjugated  remaining  elements; 

prod2  =  prodl*  prodl*  . prodl*  {2)  .  (3.3) 

The  two  sequences  are  then  concatenated  and  a  zero  is  inserted  between  them  to  form  a 
proper  correlation  function  having  a  real  valued  Fourier  transform.  The  1-1/2  Djps 
estimate  can  then  be  obtained  by  multiplying  the  real  part  of  the  transformed  sequence  by 
one-half.  However,  since  the  magnitude  squared  estimate  is  required  later  in  this  study 
to  compare  the  1-1/2  and  spectrogram  representations,  ONE_HALF  computes  the 
magnitude  vice  the  real  part  of  the  transform: 

1-1/2  “  2^  ^  \prodl  0  prod^}  i  (3  4) 

where  3  denotes  Fourier  transform. 

SPECTRO  is  the  extrinsic  Matlab  function  included  in  Appendix  C  that  implements 
the  spectrogram.  It  has  the  same  input  and  output  format  as  ONE_HALF.  The  two 
functions  are  equivalent  with  respect  to  zero  padding  and  transform  length  so  that  fair 
comparisons  can  be  made  between  the  two  methods.  At  each  step  the  function  computes 
the  periodogram: 

winlen  - 1  ^ 

xin)w(n)exp{-j(on)  ;  (3.5) 

«  •  0 

where  the  data  segment.  jc(n),  is  padded  with  zeros  to  the  transform  length  of  winlen  if 
necessary,  and  w{n)  is  an  appropriate  window  fimction.  A  rectangular  window  (i.e., 
w{n)  =  1  over  the  support  of  n)  is  assumed  in  the  estimate  variance  calculations  to  follow. 
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B.  NOISE  ONLY  PERFORMANCE  COMPARISON 


The  variance  of  the  1-1/2  Djpj  estimate  is  derived  for  a  zero  mean  white  Gaussian 


input.  The  details  of  this  calculation  are  contained  in  Appendix  A.  If  the  input  is  real, 
the  theoretical  variance  based  on  (2.24)  is  given  by. 


var  { 1-1/2  = 


5^  +  6winlen  2\^ 


(oi)  ; 


(3.6) 


where  winlen  is  the  length  of  the  data  window  and  G^is  the  variance  of  the  input.  For  a 
complex  input  the  variance  is 

var  {1-1/2  Dips)  =  [winlen  +  2]  (3.7) 

The  variance  of  the  periodogram  as  defined  by  (3.5)  with  a  rectangular  window  is  [Refs. 

2.11]: 

-y  ->2 

var  {Per}  s winlen  .  (3.8) 

Both  estimate  variances  are  dependent  upon  the  length  of  the  data  window  and  the  input 
variance.  For  a  given  input  variance  and  increasing  window  length,  the  periodogram ’s 
variance  increases  more  rapidly  than  does  the  estimate  variance  of  1-1/2  Djp^  because  the 

periodogram ’s  variance  is  proportional  winlen^  rather  than  winlen.  However,  since  the 
1-1/2  Dips  variance  is  a  function  of  the  input  variance  cubed  and  the  periodogram ’s  vari¬ 
ance  depends  upon  input  variance  squared,  1-1/2  Djps’s  variance  increases  quicker  than 
the  periodogram ’s  for  a  fixed  window  length  and  increasing  input  variance. 

Computer  simulations  confirm  the  theoretical  variance  expressions  (3.6)  and  (3.8), 
and  the  expected  trends  as  window  length  and  input  variance  are  changed.  Real  valued 
white  Gaussian  noise  is  processed  by  ONE_HALF  employing  a  rectangular  window  and  a 
step  increment  set  equal  to  the  window  length  so  that  no  overlapping  of  data  segments 
occurs.  Since  SPECTRO  is  written  for  a  fixed  step  increment  of  one,  the  periodogram  is 
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estimated  directly  from  (3.5)  in  the  Matlab  script  file  used  to  perform  the  simulations. 
Simulation  parameters,  measured  estimate  variances,  and  theoretically  expected  variances 
are  shown  in  Table  3.1.  Measured  input  variances  are  used  to  compute  the  theoretical 
variances.  There  is  good  agreement  between  theoretical  and  expected  values.  The 
methods  cannot  be  compared  number  for  number  to  each  other  because  the  periodogram 
variance  is  a  magnitude  squared  quantity  and  1*1/2  ^ips  is  a  magnitude  quantity. 
However,  comparisons  between  the  methods  can  be  made  regarding  their  relative 
responses  to  changing  window  length  and  input  variance.  For  a  window  length  of  256, 
1-1/2  Djps’s  variance  increases  by  a  factor  of  about  342  when  input  variance  is  increased 
from  one  to  seven.  For  the  same  variance  change,  the  periodogram ’s  variance  increases 
by  a  considerably  smaller  factor  of  about  49.  When  input  variance  is  three  and  window 
length  increases  from  128  to  512  the  periodogram’s  variance  increases  by  a  factor  of  15 
while  1-1/2  Dips  increases  by  a  factor  of  three. 


TABLE  3.1:  THEORETICAL  VS.  MEASURED  M/2  Dn»s  AND  PERIODOGRAM 
VARUNCES  FOR  VARIOUS  WINDOW  LENGTHS  AND  INPUT  VARIANCES 


input 

variance 

window 

length 

1-1/2  Dips 

theoretic^ 

1-1/2  Dips 

measured 

periodogram 

theoretical 

periodogram 

measured 

variance 

variance 

variance 

variance 

1 

128 

204.3 

205.5 

20,456 

19,487 

1 

256 

394.2 

377.8 

77,549 

93,217 

1 

512 

779.8 

707.6 

303,490 

297,930 

3 

128 

5,515.8 

5,549.5 

184,100 

175,380 

3 

256 

10,644 

10,201 

697,940 

838,960 

3 

512 

21,056 

19,107 

2,731,400 

2,681,300 

7 

128 

70,071 

70,499 

1,002,300 

954,850 

7 

256 

135,220 

129,590 

3,799,900 

4,567,600 

7 

512 

267,490 

242,720 

14,871,000 

14,598,000 
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The  1-1/2  Djps  and  spectrogram  time-frequency  representations  of  a  white  Gaussian 
noise  input  are  shown  in  Figure  3.1.  ONE_HALF  generates  the  1-1/2  Djpg  estimate  and 
SPECTRO  produces  the  spectrogram.  Both  representations  are  formed  by  applying  a 
rectangular  window  to  a  single  128  noise  sample  realization  using  a  64  point  window  length 
and  a  step  increment  of  one.  The  most  striking  difference  between  the  two  surfaces  is  that 
1-1/2  Dips  consists  of  several  short  time  duration  spikes  while  the  spectrogram  exhibits  a 
more  continuous  and  smoother  surface.  The  spikier  appearance  of  the  1-1/2  Djpj  surface 
indicates  that  the  1-1/2  Djpj  estimate  de-correlates  faster  than  the  spectrogram.  As  a 
consequence,  the  typical  1-1/2  Djps  frequency  bin  contains  more  independent  segments 
than  does  the  typical  spectrogram  bin.  This  is  signiAcant  because  a  larger  reduction  in 
variance  is  realizable  from  a  greater  number  of  independent  segments  when  the  power  in  a 
bin  is  summed. 

The  number  of  independent  segments  in  a  bin  is  estimated  by  applying  a  relation  used 
to  determine  the  degrees  of  freedom  for  a  chi-squared  distribution.  Given  a  collection  of 
n  normally  distributed  zero  mean  unit  variance  random  variables,  the  degrees  of  freedom 
(DF)  represent  the  number  of  independent  random  variables  contained  in  the 
collection.  Degrees  of  freedom  are  found  from  the  mean  and  variance  of  the  sum  of  the 
random  variables,  s  =  0:^  +  0:2  +  .. .  +  using  the  relation: 


where  and  are  the  mean  and  variance  respectively,  of  s.  [Ref.  7] 

Ten  spectrogram  and  1-1/2  ^ips  representations  for  different  realizations  of  real  white 
unit  variance  Gaussian  noise  are  used  to  estimate  the  degrees  of  freedom  with  (3.9).  The 
degrees  of  freedom  are  calculated  for  the  individual  time-frequency  cell  (i.e.,  the  estimate’s 
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frequency  bln 

(C)  (d) 

Figure  3.1:  Time-frequency  representations  of  unit  variance  white  Gaussian  noise; 
(a)  1-1/2  D|pg  mesh  plot,  (b)  1-1/2  contour  plot,  (c)  spectrogram  mesh  plot,  and 

(d)  spectrogram  contour  plot. 
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value  in  one  frequency  bin  at  one  time  step),  as  well  as  for  the  frequency  bin.  The 
degrees  of  freedom  for  a  cell  is  based  upon  the  mean  and  variance  of  all  the  values  in  the 
frequency  bin.  The  degrees  of  freedom  for  a  bin  is  based  upon  the  mean  and  variance  of 
the  average  frequency  bin  values  for  all  the  bins  in  the  surface.  The  final  degrees  of 
freedom  estimates  are  obtained  by  averaging  the  degi  ees  of  freedom  computed  for  each  of 
the  ten  realizations.  Table  3.2  summarizes  the  results  of  the  simuktions  for  various 
window  lengths  and  window  functions.  In  all  cases.  SPECTRO  and  ONE_HALF  use  a 
step  increment  of  one,  and  the  data  sequence  is  twice  the  size  of  the  window  length.  The 
degrees  of  freedom  for  a  cell  in  the  1-1/2  Djp^  representation  is  approximately  0.6  while 

the  degrees  of  freedom  for  the  spectrogram  cell  is  about  four.  The  degrees  of  freedom 
for  a  bin  depends  upon  the  window  function.  With  a  rectangular  window,  a  1-1/2  Djpj 
bin  has  about  17.5  degiees  of  freedom  and  a  spectrogram  bin  has  almost  six  degrees  of 
freedom.  If  a  Hamming  window  is  applied  to  the  data,  the  spectrogram  and  1-1/2  Djpj 
bin  degrees  of  freedom  figures  increase  to  about  eight  and  27,  respectively. 


TABLE  3.2:  DEGREES  OF  FREEDOM  FOR  1-1/2  Dips  AND  SPECTROGRAM 

ESTIMATES 


winlen 

window 

frmction 

ONEHALF 

DFceii 

ONEHALF 

DFbin 

SPECTRO 

DFceU 

SPECTRO 

DPbin 

64 

rectangular 

0.6582 

16.3737 

4.1109 

5.8239 

64 

Hamming 

0.6185 

26.9717 

4.4057 

7.7394 

128 

rectangular 

0.6326 

18.7693 

4.0822 

5.9613 

128 

Hamming 

0.5974 

27.5028 

4.1742 

8.0081 

The  data  sequence  length  divided  by  the  bin  degrees  of  freedom  indicate  the  size  of 
the  independent  segments  in  a  typical  frequency  bin.  For  example,  the  size  of  an 
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independent  segment  in  the  1-1/2  Djps  estimate  based  on  a  window  length  of  64,  a  128  point 

sequence,  and  a  rectangular  window  is  128/16.3737  s  8.  For  the  same  parameters,  the 
spectrogram  independent  bin  segment  size  is  about  22.  These  figures  seem  reasonable 
based  upon  “eyeball”  estimates  of  inde{>endent  bin  segment  size  in  Figure  3. 1  contour  plots. 

When  an  estimate’s  h’equency  bin  values  are  averaged,  the  anticipated  variance 
reduction  factor  is  given  by  the  ratio  of  that  estimate’s  DF^in  to  DFj^  values.  The 
expected  variance  reduction  factor  for  1-1/2  ^ips  is  about  29  for  a  rectangular  window  and 
45  for  a  Hamming  window.  The  spectrogram  has  an  expected  variance  reduction  factor 
of  1 .5  for  a  rectangular  window  and  two  for  a  Hamming  window. 

Computer  simulations  using  ONE_HALF  and  SPECTRO  confirm  that  an  exploitable 
difference  in  variance  reduction  capability  exists  between  1-1/2  Djps  and  the  spectrogram 
when  frequency  bin  averages  of  a  time-frequency  surface  are  computed.  Gaussian  white 
noise  is  simulated  to  serve  as  a  noise  only,  real  input  signal.  Both  methods  employ  a 
rectangular  window  on  the  same  input  signal.  The  step  increment  is  specified  to  be  one 
for  each  method  and  the  data  window  length  is  consistently  set  equal  to  half  of  the  data 
sequence  length.  For  the  first  set  of  simulations,  noise  variance  is  fixed  while  the 
transform  length  is  varied.  The  simulations  require  that  each  method  generate  a  time- 
frequency  representation  for  fifteen  different  input  sequences .  The  fifteen  representations 
are  then  averaged  to  form  one  time-frequency  surface.  For  this,  and  all  remaining 
simulations  performed  in  this  study,  the  averaged  1-1/2  Djpj  surface  is  squared  so  that  its 
amplitude  can  be  fairly  compared  to  the  spectrogram’s  amplitude.  As  a  consequence  of 
computer  program  design,  data  sequence  length,  and  choice  of  step  size  the  first  and  last 
winlenP.  time  steps  are  calculated  using  zero-padded  data  segments.  These  time  steps  are 
deleted  and  the  average  value  of  each  frequency  bin  on  the  average  surface  is  calculated. 
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Figures  3.2  through  3.6  show  these  average  frequency  bin  values  for  an  input  noise 
variance  of  one  and  window  lengths  of  32,  64,  128,  256,  and  512,  respectively.  The  bin 
corresponding  to  dc,  and  the  bin  corresponding  to  the  fold-over  frequency  are  omitted  from 
the  plots.  To  allow  for  a  meaningful  comparison,  the  mean  frequency  bin  average  is 
subtracted  from  each  plot.  Examination  of  Figures  3.2  through  3.6  indicate  that  1-1/2  Djpj 
exhibits  a  lesser  output  variation  than  does  the  spectrogram  as  the  transform  length 
increases.  This  observation  is  confirmed  in  Figure  3.7  which  plots  the  measured  output 
variance  versus  transform  length  for  the  first  simulation  set. 


1>1/2  Dips 


Spsct  rog  ram 


Figure  3.2:  Averaged  output  variance  for  winlen  =  32  and  =  1 . 
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Figure  3.4:  Averaged  output  variance  for  winlen  =  128  and 
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Figure  3.7:  Spectrogram  (dashed)  and  1-1/2  Djpg  (solid),  output  variance  vs.  window  length 

{oral  =  1. 


The  simulation  was  repeated  for  different  input  noise  variances  and  a  fixed  window 
length.  Figure  3.8  shows  the  relationship  between  output  variance  and  input  variance  for 
a  window  length  of  128.  The  point  at  which  the  two  curves  in  Figure  3.8  intersect  can  be 
considered  a  performance  breakpoint.  At  input  variances  below  the  breakpoint  value  (i.e. , 
about  1.5  from  Figure  3.8),  the  1-1/2  Djpj  estimate  has  less  variation  than  the 
spectrogram.  When  the  input  variance  is  larger  than  the  breakpoint  value  the  spectrogram 
exhibits  less  variation  than  1-1/2  Djpg.  Other  simulations  show  that  breakpoints  for 
window  lengths  ranging  from  32  to  512,  although  gradually  increasing  as  window  length 
increases,  all  correspond  to  an  input  noise  variance  of  approximately  1.5. 

Figure  3.8  shows  that  the  spectrogram  has  a  variance  of  about  1000  when  the  input 
variance  is  about  1.3.  Considering  the  processing  scheme,  the  expected  variance  is  the 
theoretical  estimate  variance  given  by  (3.8)  reduced  by  two  factors.  The  first  factor  is 
equal  to  the  number  of  time-frequency  surfaces  averaged  to  obtain  a  representative  surface. 
The  second  factor  is  the  degrees  of  freedom  based  variance  reduction  factor  estimated 
earlier.  Therefore,  the  expected  spectrogram  output  variance  is  given  by: 
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7  7  2 

winlen  (a^ 

VT  {Specro}  u 

__  128^  •  1.32 

15  •  1.5 
=  1231; 

where  the  numerator  is  (3.8),  nsurf  is  the  number  of  time-frequency  surfaces  averaged  in 
the  simulation,  and  VRF  is  the  degrees  of  freedom  variance  reduction  factor  for  a 
rectangular  window.  The  expected  and  actual  variances  for  the  spectrogram  agree  fairly 
well  considering  that  the  comparison  is  based  upon  approximate  values. 

The  same  comparison  is  made  for  1-1/2  Djp^.  According  to  Figure  3.8,  an  input 
variance  of  approximately  1.45  results  in  a  1-1/2  Djpg  estimate  variance  of  about 
1000.  The  theoretical  variance  given  by  (3.6)  is  reduced  by  the  same  factor  of  15  used  in 
the  spectrogram  calculation  since  this  factor  depends  only  upon  the  number  of  time- 
frequency  surfaces  averaged.  The  theoretical  variance  is  also  reduced  by  a  degrees  of 
freedom  variance  reduction  factor  that  is  larger  than  the  spectrogram’s  (29  vice  1.5).  To 
accoimt  for  squaring  the  output  of  ONE.HALF  in  the  simulations,  (3.6)  is  also 
squared.  The  expected  1-1/2  Djps  variance  based  on  the  processing  scheme  is: 


var  { 1-1/2  Dips)  s 


[54  +  6winlen]^{ol) 

M _ •* 

4^  •  nsurf  -  VRF 
[54 6- 128]  2  (1.45)' 
4^  •  15  •  29 

=  902. 


The  expected  variance  agrees  fairly  well  with  the  simulation  variance  for  the  1-1/2  Dips 
estimate  considering  that  approximate  values  are  being  compared. 
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Figure  3.8;  Spectrogram  (dashed)  and  1-1/2  Djpg  (solid),  output  variance  vs.  input  variance 

for  winlen  =  128. 


C.  SIGNAL  ONLY  PERFORMANCE  COMPARISON 

Appendix  B  details  the  calculations  made  using  (2.24)  to  determine  the  theoretical 
1-1/2  Djpj  estimate  for  a  signal  only.  For  a  real  valued  sinusoid  given  by 

x(n)  s  cos(2n^/i),  (3.10) 

the  theoretical  estimate  is: 

A  ^  m 

1-1/2  Djp5(Ai,©)  =  2^08  (2n^/i)  {6(a)-m) +6[a)- (N-m)] } .  (3.11) 

The  periodogram  is  defined  in  (3.5).  For  the  same  signal,  the  theoretical  value  of  the 
periodogram  is  the  magnitude  squared  of  the  Fourier  transform  of  the  input  signal: 

N  ^ 

Perin,(o)  =  (2)  {6(©-m) +6[©- (N-m)]  }  (3.12) 

If  the  signal  consists  of  a  single  complex  sinusoid, 

x(n)  =  expijln^n),  (3.13) 
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the  theoretical  1*1/2  Dips  estimate  is  given  by; 

1-1/2  Djpj (/I,©)  =  Ncos(2n^/i)6(©--m) ;  (3.14) 

and  the  theoretical  value  of  the  periodogram  is: 

Per(/i,  ©)  =  (3.15) 

Results  of  computer  simulations  using  SPECTRO  and  ONE_HALF  validate  the 
theoretical  calculations  for  the  periodogram  and  1-1/2  Dips*  respectively.  For  the 
simulations,  a  128  point  data  window  is  stepped  through  a  256  point  data  sequence  one 
point  at  a  time.  Both  methods  apply  a  rectangular  window  function  to  the  data.  The  first 
and  last  winlenP.  time  steps  are  discarded  before  computing  any  averages  or  plotting  any 
results.  The  1-1/2  Djp^  estimate  is  squared  to  allow  a  fair  comparison  with  the 
spectrogram.  Figure  3.9  shows  the  mesh  plot  and  corresponding  contour  plot  of  the 
1-1/2  Djps  estimate  for  the  real  signal.  The  sinusoidal  signal  is  centered  in  frequency  bin 
20.  Bin  20  and  the  average  of  each  frequency  bin  are  plotted  in  Figure  3.10(a)  and 
Figure  3.10(b),  respectively.  The  maximum  power  at  any  one  time  step  is  approximately 
4042.  This  agrees  fairly  well  with  the  magnitude  squared  value  of  (3.11)  which  is  equal 
to  half  of  the  window  length  squared,  or  4096.  The  intra-ridge  modulation  exhibited  by 
bin  20  in  Figures  3.9  and  3.10(a)  results  from  the  cosine  cubed  modulating  term  in  (3.11). 
The  average  power  of  bin  20  is  approximately  1280  as  shown  in  Figure  3. 10(b). 

Figures  3.11  and  3.12  use  the  same  arrangement  of  plots  to  display  spectrogram  results 
for  the  real  signal.  The  maximum  power  agrees  exactly  with  the  theoretical  value  of 4096 
calculated  using  (3.12).  Figure  3.12  shows  that  the  average  power  of  bin  20  is  equal  to 
the  maximum  power.  The  maximum  power  of  each  method  is  about  the  same,  however, 
the  average  power  of  1-1/2  Djpj  is  approximately  one- third  as  much  as  the  spectrogram’s 
average  power  due  to  the  intra-ridge  modulation  effect. 
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Figure  3.12:  Spectrogram  signal  power  in  bin  20  (real  sinusoid). 


When  the  same  simulations  are  run  with  a  complex  sinusoidal  input,  ONE_HALF 
produces  a  maximum  signal  power  ^^proximately  equal  to  the  length  of  the  data  window 
squared  which  agrees  with  the  squared  magnitude  of  (3.14).  As  a  result,  ONE_HALF’s 
maximum  power  is  about  equal  to  the  spectrogram’s  for  both  real  valued  and  complex 
valued  signals.  However,  for  the  complex  valued  sinusoid  the  average  power  in  bin  20  of 
the  1-1/2  D|ps  estimate  is  about  half  that  of  the  spectrogram’s  rather  than  a  third  as  it  was 
for  the  real  sinusoid.  This  is  a  consequence  of  the  different  modulating  terms  in  (3. 1 1)  and 
(3.14).  The  real  signal  is  modulated  by  a  cosine  cubed  function  which  has  a  magnitude 
squared  average  value  of  about  0.345  for  a  unit  amplitude  cosine  over  one  period.  In 
contrast,  the  complex  valued  signal  is  modulated  by  a  cosine  function  having  a  magnitude 
squared  average  value  of  about  0.52. 
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techniques  forms  a  figure  of  merit  used  to  compare  methods.  To  compare  1-1/2  and 
the  spectrogram  using  (3.16)  the  processing  gain  ratio  is  defined  as: 


PG 


SNR  - 


f*^^^ONEHALF 

PSNR^p^qjhq 


When  (3.17)  is  the  measure  of  interest,  the  processing  gain  ratio  is  given  by: 


PG 


DIF  - 


^^^^ONEHALF 

PDIF^pectro 


(3.18) 


(3.19) 


Computer  simulations  are  used  to  measure  the  processing  gain  for  various  window 
lengths  and  input  noise  variances.  The  input  consists  of  a  fixed  amplitude  real  sinusoid 
centered  in  frequency  bin  20  and  mixed  with  Gaussian  white  noise.  For  this  scenario  the 
input  signal-to-noise  ratio  (SNR)  is  given  by: 


SNR  =  lOlog 


(3.20) 


Fifteen  realizations  of  each  simulation  are  averaged  in  order  to  obtain  representative 
results.  As  before,  the  input  data  sequence  is  twice  the  length  of  the  data  window,  a 
rectangular  window  is  applied  to  the  data,  and  the  data  window  is  stepped  through  the 
input  sequence  one  point  at  a  time.  Figure  3.13  shows  the  resulting  plots  for  a  data 
window  length  of  128,  and  an  input  SNR  of  -3dB.  The  fifteen-realization  averaged 
1-1/2  Dips  time-ffequency  surface  is  shown  in  Figure  3.13(a)  and  its  corresponding 
frequency  bin  average  is  displayed  in  Figure  3.13(b).  Figures  3. 13(c)  and  3. 13(d)  are  the 
averaged  spectrogram  time-ffequency  surface,  and  its  frequency  bin  averages, 
respectively.  Figures  3.14,  3.15,  3.16,  3.17,  and  3.18  display  the  results  for  the  same 
window  length  but  for  different  input  SNRs  of  about  -lOdB,  -15dB,  -17dB,  -18.5  dB,  and 
-19.5  dB,  respectively. 
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Processing  gain  comparisons  using  (3.18)  and  (3.19)  were  made  of  the  frequency  bin 
averages  as  a  function  of  SNR  for  the  signals  shown  in  Figures  3.13  through  3.18.  The 
average  value  of  the  signal  bin  was  used  for  P^ig  to  compute  PSNR  and  PDIF.  The  mean, 
standard  deviation,  and  variance  of  the  noise  were  computed  using  all  bins  except  the  one 
containing  the  signal.  Table  3.3  summarizes  the  results  of  these  processing  gain 
computations. 

The  PGsiiqi  performance  measure  is  consistent  with  results  obtained  from  the  noise 
only  simulations  and  the  signal  only  simulations.  As  expected  from  these  earlier 
simulations,  PG^nr  is  greater  than  one  (indicating  that  1-1/2  ^ips  performs  better  than  the 
spectrogram)  when  the  input  variance  is  smaller  than  about  1.5.  This  is  agrees  with  the 
breakpoint  determined  from  the  analysis  of  Figure  3.8  earlier.  Based  upon  the  time- 
frequency  representations  however,  PGdip  appears  to  be  a  more  useful  indicator  of 
detection  performance  over  a  much  wider  range  of  input  SNRs. 

The  PGdip  column  of  Table  3 .3  shows  that  the  1  - 1/2  Djpj  frequency  bin  average  signal 
peak  estimate,  when  measured  in  terms  of  noise  standard  deviations  above  the  mean  noise 
level,  is  about  1.6  times  higher  than  it  is  for  the  spectrogram  estimate  when  the  input  SNR 
is  about -3dB.  As  input  SNR  decreases  the  relative  processing  gain  also  decreases.  Both 
methods  are  equivalent  when  input  SNR  is  about  -19dB.  The  spectrogram  becomes  the 
better  estimate  by  this  standard  at  SNRs  lower  than  -19dB  for  this  window  length. 
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Figure  3.14:  Sinusoid  in  bin  20  with  SNR  s  -lOdB;  (a)  1-1/2  Djpj  time-frequency 
representation,  (b)  1-1/2  Djpj  bin  averages,  (c)  spectrogram  time-frequency  representation, 

and  (d)  spectrogram  bin  averages. 


Figure  3.15:  Sinusoid  in  bin  20  with  SNR  s  -15dB;  (a)  1-1/2  Djps  time-frequency 
representation,  (b)  1-1/2  Djpj  bin  averages,  (c)  spectrogram  time-frequency  representation, 

and  (d)  spectrogram  bin  averages. 
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Figure  3.16;  Sinusoid  in  bin  20  with  SNR  s  -17dB;  (a)  1-1/2  Djpj  time-frequency 
representation,  (b)  1-1/2  averages,  (c)  spectrogram  time-frequency  representation, 

and  (d)  spectrogram  bin  averages. 


\5 


tr«qu«ncy  bln 


frequency  bln 


(c)  (d) 

Figure  3.17:  Sinusoid  in  bin  20  with  SNR  s  -18.5dB;  (a)  1-1/2  time-frequency 
representation,  (b)  1-1/2  bin  averages,  (c)  spectrogram  time-frequency  representation, 

and  (d)  spectrogram  bin  averages. 
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Figure  3.18:  Sinusoid  in  bin  20  with  SNR  ss  -19.5dB;  (a)  1-1/2  Djpj  time-frequency 
representation,  (b)  1-1/2  bin  averages,  (c)  spectrogram  time-frequency  representation, 

and  (d)  spectrogram  bin  averages. 


TABLE  3.3:  RELATIVE  PROCESSING  GAINS  FOR  VARYING  INPUT  SNR  AND 
A  FIXED  DATA  WINDOW  LENGTH  OF  128 


input  SNR 
(-dB) 

FGsnr 

FDI^onehalf 

I^l^SPECTRO 

conespooding 
figure  no. 

3.01 

1.4270 

217.8826 

132.4518 

1.6450 

3.12 

10.00 

<  1 

37.9275 

29.3067 

1.2942 

3.13 

14.77 

<1 

12.8635 

11.2700 

1.1414 

3.14 

16.99 

<  1 

8.0222 

7.4846 

1.0718 

3.15 

18.45 

<  1 

5.9961 

5.8097 

1.0269 

3.16 

19.54 

<  1 

4.8344 

4.8552 

0.9957 

3.17 

Results  of  signal-in-noise  simulations  run  for  a  smaller  64  point  data  window  and  a 
larger  256  point  data  window  are  consistent  with  those  obtained  for  the  128  point 
window.  However,  the  range  of  input  SNRs  over  which  the  1-1/2  Djps  method  yields  a 
larger  processing  gain  than  does  the  spectrogram  depends  upon  window  length.  For  the 
64  point  window,  1-1/2  D|ps  yields  a  larger  processing  gain  than  the  spectrogt  un  only  for 
input  SNRs  greater  than  about  - 1 2 .5dB .  As  already  demonstrated,  for  a  1 28  point  window 
1-1/2  Dips  ^  larger  processing  gain  than  the  spectrogram  for  input  SNRs  greater  than 
about  -19dB.  Simulations  using  a  256  point  window  showed  that  1-1/2  Djps  performs 
better  than  the  spectrogram  in  terms  of  processing  gain  even  when  the  input  SNR  is  less 
than  -21dB  and  the  signal  cannot  be  readily  discerned  from  the  noise  in  the  averaged 
frequency  bin  plots.  In  effect,  for  this  type  of  signal  and  processing  scheme  using  a  256 
point  data  window,  1- 1/2  Dip^’s  processing  gain  is  between  1.15  and  1 .64  times  larger  than 
the  spectrogram’s  for  input  SNRs  between  -21dB  and  -3dB.  The  simulations  also  indicate 
that  the  processing  gain  performance  of  1-1/2  Djp^  improves  relative  to  the  spectrogram 
with  increasing  window  length. 
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IV.  SIMULATION  RESULTS  AND  ANALYSIS 


Simulation  results  are  presented  to  show  how  1-1/2  Dips.  the  bispectrum  methods. 
IHOMS,  and  the  spectrogram  deal  with  three  types  of  signals.  The  spectrogram, 
bispectnun,  and  1-1/2  Djpj  methods  arc  first  compared  using  a  signal  consisting  of  multiple 
unrelated  stationary  sinusoids  that  are  mixed  with  white  Gaussian  noise.  Next,  the 
bispectral  representation  of  noise  free  signals  containing  harmonically  related  stationary 
components  is  studied.  Lastly,  a  received  signal  containing  multiple  linear  FM 
components  mixed  with  white  Gaussian  noise  is  processed  using  IHOMS,  1-1/2  Djps,  and 
the  spectrogram. 

A.  STATIONARY  SINUSOIDS 

The  spectrogram  is  computed  using  the  Matlab  function  SPECTRO,  and  the  1- 1/2  Djpj 
method  is  implemented  via  ONE_HALF.  Their  time-frequency  representations  for  a 
single  sinusoidal  signal  are  discussed  in  Chapter  m. 

The  bispectrum  methods  produce  a  frequency  frequency  representation  instead  of  a 
time-frequency  representation.  The  direct  method  for  computing  the  bispectrum  is 
realized  through  the  Hi-Spec  Matlab  function  BISPEC_D  of  [Ref.  14].  This  function 
returns  the  complex  bispectrum  matrix  and  a  frequency  axis  labeling  vector.  The  user 
specifies  an  input  data  sequence,  the  desired  Fast  Fourier  Transform  (FFT)  size,  the  number 
of  samples  per  segment,  the  amount  of  overlap  between  data  segments,  and  the  desired 
form  of  the  frequency  smoothing  window.  The  extrinsic  Matlab  fxmction  INDBIS 
contained  in  Appendix  C  implements  the  indirect  method  of  calculating  the  bispectrum.  It 
returns  the  complex  bispectrum  of  the  signal  supplied  by  the  user.  The  FFT  size,  the 
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number  of  data  samples  per  segment,  and  the  desired  tri-correlation  window  function  are 
specified  by  the  user.  The  three  window  options  available  are  the  unit  hexagonal  window, 
the  Parzen  window,  and  the  optimum  window  (also  known  as  the  minimum  bias  supremum 
window). 

Both  bispectrum  methods  are  used  to  generate  bispectral  representations  of  a  noise 
free  sinusoid  located  at  frequency  bin  20.  The  input  sequence  consists  of  256 
samples.  Both  methods  apply  a  128  point  FFT  to  each  of  four  non-overlapping  data 
segments  containing  64  data  points  each.  The  indirect  method  representations  are  formed 
fi-om  tri-correlation  functions  based  upon  31  lagr.  Figure  4.1  shows  the  full  bispectrum 
plane  representation  of  a  real  sinusoid  generated  by  the  direct  method  without  a 
window  The  sinusoid  is  indicated  by  the  peak  in  the  first  quadrant  located  at 
(k^,k2)  =  (20,  20)  .and  associated  symmetry  peaks  like  the  one  in  the  third  quadrant  at 
iki,k2)  *  (-20,-20). 

The  representation  of  a  real  sinusoid  at  bin  20,  formed  by  the  indirect  method  using  a 
unit  hexagonal  window  is  displayed  in  Figure  4.2.  The  rectangular  window  used  by  both 
ONE_HALF  and  SPECTRO,  the  unit  hexagonal  window  used  by  INDBIS,  and  the 
unwindowed  BISPEC_D  are  essentially  equivalent  window  functions  so  that 
representations  formed  using  them  are  considered  unwindowed.  Two  closely  spaced 
peaks  are  discernible  at  the  signal  locations  in  the  indirect  method’s  representation.  This 
is  a  characteristic  of  this  method  which  depends  upon  the  number  of  lags  used  to  compute 
the  tri-coirelation  fimction.  The  two  peaks  at  each  location  become  less  distinct  as  the 
number  of  lags  used  to  compute  the  tri-correlation  sequence  is  increased.  Another 
difference  between  Figure  4.1  and  4.2  is  that  the  indirect  representation  exhibits  peaks  at 
(ki,k2)  -  (20,-20)  and  (^|,it2)  =  (-20,20)  which  are  not  noticeable  in  the  direct 
method’s  representation.  It  turns  out  that  both  of  the  methods  are  sensitive  to  the  phase 
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of  a  real  signal,  and  that  the  methods  produce  slightly  different  representations  for  the  same 
phase.  The  signal  used  to  create  the  bispectrum  in  Figure  4.1  has  a  different  phase  than 
the  signal  which  the  bispectrum  in  Figure  4.2  is  based  upon.  Both  phase  values  are 
specially  chosen  to  chosen  to  minimize  any  peaks  on  the  it  j  =  0  and  kj  =  0  axes.  In  the 
case  of  the  direct  method’s  representation,  this  also  minimized  the  peaks  at 
(jt^,  k2)  -  (20,  -20)  and  (itj,  k2)  =  (-20, 20) .  Bispectral  peaks  located  on  either 
the  ki  or  the  it2  axis  in  subsequent  plots  are  referred  to  as  zero  axis  peaks. 

The  unwindowed  full  bispectrum  representations  for  a  complex  sinusoid  obtained  by 
the  direct  method  and  indirect  method  are  shown  in  Figure  4.3  and  4.4,  respectively.  Like 
the  real  signal  representation,  the  complex  signal  produces  a  peak  at  (k^,  ^2)  =  (20, 20) . 
The  complex  signal’s  bispectrum  does  not  exhibit  as  many  symmetry  peaks  as  does  the  real 
signal’s  representation  because  the  bispectrum  of  a  complex  signal  is  symmetric  about  only 
one  symmetry  line  and  in  only  one  quadrant.  Which  line,  and  which  quadrant  are 
determined  by  the  conjugation  scheme  used  to  compute  the  bispectrum  [Ref.  9].  For  the 
equivalent  schemes  used  by  INDBIS  and  BISPEC_D  which  are  described  in  Chapter  H, 
symmetry  exists  about  the  k^  =  k2  line  in  the  first  quadrant.  The  peaks  at 
{ki,k2)  =  (0,20)  and  ik^,k2)  =  (20,0)  in  both  representations  are  zero  axis 
peaks.  Unlike  real  signal  representations,  zero  axis  peaks  in  the  bispectrum  of  a  complex 
signal  have  a  consistent  amplitude  regardless  of  signal  phase. 

The  Brst  bispectral  quadrant  is  essentially  the  same  for  either  a  real  or  a  complex 
sinusoid,  and  is  sufficient  for  measuring  a  sinusoidal  signal’s  frequency.  The  only 
additional  information  not  contained  in  the  first  quadrant  is  whether  or  not  the  signal  is 
complex  or  real.  The  remaining  bispectrum  figures  in  this  thesis  will  display  just  the  first 
quadrant. 
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The  next  four  examples  of  single  real  valued  sinusoid  bispectrums  are  intended  to 
show  what  windowing  can  do  for  each  method.  At  the  same  time,  the  effects  of  signal 
phase  on  the  bispectral  rq)resentations  are  noted.  Figure  4.S  shows  the  bispectrum  as 

computed  by  the  direct  method  and  smoothed  by  a  Rao-Gabr  window  over  5  adjacent 
frequency  points.  The  signal  is  a  real  sinusoid  located  at  bin  20  with  a  phase  angle  of  0.484 
radians.  The  same  processing  scheme  yields  the  representation  shown  in  Figure  4.6  if  the 
signal  phase  is  changed  to  4.3151  radians.  The  indirect  method,  employing  a  Parzen 
window  on  a  real  signal  with  a  phase  equal  to  0.484  radians  produces  the  representation 
shown  in  Figure  4.7.  When  the  signal  with  a  phase  of  4.3151  radians  is  processed  by  the 
indirect  method  using  an  optimum  window  (the  optimum  window  and  Parzen  window  are 
not  significantly  different  for  this  type  of  signal),  it  produces  the  bispectrum  shown  in 
Figure  4.8.  Windowing  helped  to  make  the  two  distinct  peaks  at  ^2)  =  (20, 20) 
appear  as  one  in  the  indirect  representation.  The  zero  axis  peaks  of  the  indirect 
representation  incorrectly  appear  to  be  offset  from  the  zero  axes  because  the  fu^t  quadrant 
plots  do  not  show  the  adjacent  peaks  in  the  second  and  fourth  bispectral  quadrants.  The 
signal  frequency  is  more  accurately  and  more  clearly  depicted  by  the  direct  method’s  peak 
locations .  In  addition,  the  direct  method  is  faster  than  the  indirect  method  as  implemented 
in  this  study.  However,  if  the  signal  is  known  to  be  either  complex  or  real,  the  indirect 
method  can  be  implemented  in  such  a  way  as  to  take  advantage  of  the  bispectrum  symmetry 
in  order  to  reduce  computational  costs.  This  is  done  in  [Ref.  14]  to  implement  the  indirect 
method  for  real  signals  in  Hi-Spec’s  Matlab  function  BISPEC.I. 

The  bispectrum,  spectrogram,  and  1-1/2  Djps  methods  are  now  used  to  process  a  signal 
consisting  of  three  complex  valued  sinusoids  in  additive  white  Gaussian  noise.  The  three 
frequencies  of  the  signal  are  at  bins  10.65, 28.05,  and  54.7.  Noise  variance  is  fixed  while 
signal  amplitude  is  varied  to  achieve  different  input  SNRs.  The  phase  angle  of  each  signal 
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Figure  4.8:  Optimum  windowed  indirect  bispectrum  of  real  sinusoid  (4.3151  radians) 


component  is  distributed  uniformly  on  the  interval  (0, 2n) .  The  data  sequence  has  256 
samples.  The  1-1/2  Djpj  spectrum  is  computed  by  ONE_HALF.  the  spectrogram  is 
generated  by  SPECTRO,  and  the  bispectrum  is  calculated  by  the  direct  method  using 
BISPEC_D.  A  step  increment  of  one  is  employed  for  both  SPECTRO  and  ONE_HALF, 
while  a  50%  overlap  is  specified  for  BISPEC_D.  To  allow  for  a  fair  comparison,  all 
three  methods  are  unwindowed.  Fifteen  representations  for  each  method  are  computed 
using  different  noise  realizations.  These  surfaces  are  then  averaged  to  obtain  either  a 
representative  time-frequency  surface,  or  a  representative  frequency-frequency  surface,  as 
appropriate.  The  averaged  surface,  along  with  its  frequency  bin  average  is  displayed  for 
various  input  SNRs  and  window  length.  The  bispectral  ri^equency  bin  average  is  the 
average  power  in  each  ki  frequency  bin.  A  horizontal  line  on  the  frequency  bin  plot  is 
placed  three  noise  standard  deviations  above  the  noise  mean  to  provide  a  reference 
point.  The  noise  statistics  are  computed  after  removing  outliers  from  the  frequency  bin 
average  values.  Figure  4.9  shows  the  spectrogram  representation  for  a  128  point  window 
length,  and  an  input  SNR  of  -6dB.  The  1-1/2  Djpj  spectrum,  and  the  bispectrum  for  the 
same  input  SNR  and  window  length  are  displayed  in  Figures  4.10  and  4.11,  respectively. 
The  signal  components  are  clearly  detected  by  all  three  methods  but  the  symmetry  of  the 
bispectrum  causes  extra  peaks  that  extend  above  the  three  standard  deviation  reference 
line  in  Figure  4.11(c). 

Figures  4.12,  4.13,  and  4.14  show  the  spectrogram,  the  1-1/2  Djps  spectrum  and  the 
bispectrum,  respectively,  for  an  input  SNR  of -18dB  and  a  window  length  of  128.  At  this 
lower  SNR  it  is  difficult  to  distinguish  the  signal  components  in  the  methods’  mesh  plots. 
Figure  4.12(a),  4.13(a),  and  4.14(a).  However,  the  sinusoids  are  discernible  in  the  contour 
plots  and  bin  average  plots  of  all  three  methods  without  any  spurious  peaks  exceeding  the 
three  standard  deviation  reference  line.  The  slight  processing  gain  advantage  of  the 
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Figure  4.1 1:  Bispectrum  of  three  complex  sinusoids  at  an  input  SNR  of  •6dB.  The 
representation  is  formed  using  a  128  point  data  window. 
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1-1/2  Djpj  method  over  the  spectrogram  at  this  SNR  and  window  length  is  barely 
noticeable  in  Figures  4.12(c)  and  4.13(c).  But,  the  signal  peaks  and  the  larger  noise 
peaks  in  the  spectrogram  plot  are  closer  to  the  three  standard  deviation  reference  line  than 
they  are  in  the  1-1/2  Djpg  plot.  Relative  to  the  spectrogram,  the  bispectrum  signal  peaks 
have  the  same  height  above  the  reference  line  but  the  noise  peaks  are  closer  to  the  line. 

Figure  4.15  displays  just  the  frequency  bin  plots  of  the  three  methods  for  an  input  SNR 
of  -21dB.  The  spectrogram  has  a  slight  advantage  over  1-1/2  Djp^  at  this  SNR  in  terms  of 
the  proximity  of  signal  peaks  and  noise  peaks  to  the  three  standard  deviation  line.  This  is 
consistent  with  Chapter  m  simulation  results  which  indicate  that  the  relative  processing 
gain  favors  the  spectrogram  below  -19.5dB  for  this  window  length.  The  bispectrum  has 
relatively  weaker  signal  peaks  and  larger  noise  peaks  than  either  of  the  other  two  methods. 

Figure  4.16  shows  the  frequency  bin  average  plots  for  the  same  SNR  of  -21dB  but  a 
larger  data  window  of  256  points.  The  1-1/2  Djpj  spectrum,  as  expected  for  this  larger 
window  length,  appears  to  perform  slightly  better  than  the  spectrogram.  However,  one 
noise  peak  at  bin  60  in  the  1-1/2  Djp^  plot  just  reaches  the  three  standard  deviation  line  in 
Figure  4.16(b).  The  bispectrum  has  strong  signal  peaks  and  relatively  low  noise  peaks 
except  for  the  one  at  about  bin  1 18  that  exceeds  the  reference  level. 

Figure  4. 17  displays  the  frequency  bin  plots  for  an  SNR  of  -26dB  and  a  window  length 
of  256.  For  all  three  methods,  only  the  first  signal  component  exceeds  the  reference  level. 
The  1-1/2  Dips  representation’s  signal  peak  has  a  slightly  larger  relative  height  above  the 
reference  line  than  either  of  the  other  methods.  No  one  method  distinguisnes  itself  with 
respect  to  suppressing  the  noise  peaks  in  this  case. 
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Figure  4.12;  Spectrogram  of  three  complex  sinusoids  at  an  input  SNR  of  -18dB.  The 
representation  is  formed  using  a  128  point  data  window. 
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Figure  4.13;  1-1/2  spectrum  of  three  complex  sinusoids  at  an  input  SNR  of  -18dB. 
The  representation  is  formed  using  a  128  point  data  window. 
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Figure  4. 14:  Bispectnun  of  three  complex  sinusoids  at  an  input  SNR  of  -18dB.  The 
representation  is  formed  using  a  128  point  data  window. 
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Figure  4.15:  Frequency  bin  averages  for  an  input  SNR  of  -21dB  and  a  128  point  window 
length;  (a)  spectrogram,  (b)  1-1/2  Djps,  and  (c)  bispectrum. 
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Figure  4.16:  Frequency  bin  averages  for  an  input  SNR  of  -216B  and  a  256  point  window 
length;  (a)  spectrogram,  (b)  1-1/2  Djpj,  and  (c)  bispectrum. 
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Figure  4.17:  Frequency  bin  averages  for  an  input  SNR  of  -26dB  and  a  256  point  window 
length;  (a)  spectrogram,  (b)  1-1/2  Djpj,  and  (c)  bispectrum. 
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B.  HARMONICALLY  RELATED  SINUSOmS 


One  of  the  main  advantages  of  higher  than  second  wder  spectra  lies  in  their  ability  to 
preserve  true  phase  information  [Ref.  1].  A  common  application  found  in  the  literature  to 
test  this  ability  is  the  detection  of  quadratic  phase  coupling  [Refs.  1,3,4,11].  Three 
sinusoids  are  quadratically  phase  coupled  when  they  are  both  harmonically  related  and 
phase  related.  Harmonically  related  means  that  one  of  the  sinusoids  has  a  frequency  equal 
to  the  sum  of  the  other  two.  Similarly,  phase  related  means  that  one  of  the  phases  is  equal 
to  the  sum  of  the  other  two.  [Ref.  1] 

Simulations  involving  BISPEC_D  show  that  the  bispectrum  is  sensitive  to  signals 
containing  harmonic  components  in  general  and  not  just  when  quadratic  phase  coupling 
occurs.  A  signal  consisting  of  two  sinusoids  is  considered  initially.  The  Hrst  sinusoid 
has  a  frequency  corresponding  to  bin  IS,  the  second  sinusoid’s  frequency  corresponds  to 
bin  40.  Both  sinusoids  have  uniform  random  phases  independently  distributed  over  the 

interval  (0, 2n) .  The  bispectrum  representation  is  generated  using  BISPEC_D  with  a  3^ 
Rao-Gabr  window,  a  64  point  data  window,  and  128  point  FFTs.  Figure  4.18  shows  the 
bispectrum  representation  of  this  signal.  Both  components  of  the  signal  are  indicated  by 
the  vertically  aligned  pair  of  peaks  in  bins  15  and  40,  and/or  by  the  horizontally  aligned 
peaks  in  k2  bins  IS  and  40.  Figure  4. 19  displa ys  the  results  for  the  same  simulation  when 
the  second  sinusoid’s  frequency  is  changed  so  that  it  now  corresponds  to  bin  30  and 
becomes  a  harmonic  of  the  first  component.  Only  the  one  peak  at  it  j  =  k2  =  IS  actually 
stands  out  on  the  bispectrum  surface.  The  bispectrum  is  blind  to  the  sinusoid  at  bin  30 
because  it  is  the  highest  harmonic  component  of  the  first  sinusoid  contained  in  the  signal. 
A  signal  comprised  of  three  harmonic  components  located  at  bins  IS,  30,  and  4S  is  now 
processed  under  two  phase  related  conditions.  For  the  first  condition,  all  the  phases  are 
i.i.d.  uniformly  over  (0, 2n) .  The  bispectrum  representation  for  this  signal  is  shown  in 
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Figure  4.20.  The  highest  harmomc  component  is  again  not  represented  by  a  significant 
peak  located  in  bin  45.  However,  the  bispectrum  does  indicate  by  its  appearance  that  the 
signal  contains  harmonic  components  based  on  the  fundamental  frequency  corresponding 
to  bin  15.  If  the  components  were  not  harmonics  of  this  fundamental  frequency  the 
bispectrum  surface  would  have  an  i^pearance  similar  to  that  previously  shown  in  Figxire 
4.11. 


(a)  (b) 

Figure  4.20:  Three  harmonic  sinusoids  in  bins  15, 30,  and  45  with  unrelated  phases. 

Figure  4.2 1  displays  the  bispectrum  for  the  three  harmonic  component  signal  under  the 
second  phase  condition.  Here  the  third  sinusoid  has  a  phase  equal  to  the  sum  of  the  other 
two  component  phases.  This  situation  satisfies  the  definition  of  quadratic  phase 
coupling.  There  is  essentially  no  difference  between  this  bispectrum  and  the  one 
generated  for  the  unrelated  phase  condition.  When  the  signal  components  are  harmonics 
of  the  same  fundamental  frequency,  the  bispectrum  representation  does  not  discern 
quadratic  phase  coupling. 
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(a)  (b) 

Figure  4.21:  Three  quadratically  phase  coupled  harmonic  sinusoids  in  bins  IS,  30,  and  45. 


If  the  signal  components  are  harmonically  related  but  are  not  harmonics  of  the  same 
fundamental  frequency,  the  bispectrum  does  distinguish  quadratic  phase  coupling.  This 
frequently  considered  situation  in  the  bispectrum  literature  is  depicted  in  Figure  4.22.  The 
signal  components  are  quadradcaUy  coupled  but  the  signal  frequencies  are  15, 25,  and  40 
rather  than  harmonics  (i.e.,  15, 30,  and  45).  This  representation  has  just  two  peaks  located 
at  the  intersection  of  the  ky  and  k2  bins  where  the  first  two  sinusoids  reside. 


(a)  (b) 

Figure  4.22:  Quadratically  phase  coupled  sinusoids  located  in  bins  15, 25,  and  45. 
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C.  MULTI-COMPONENT  LINEAR  CHIRPS 


A  computer  simulation  based  upon  Example  1  of  [Ref.  6]  is  used  to  compare  the 
processing  ability  of  IHOMS,  1-1/2  Djps,  and  the  spectrogram  when  the  received  signal 
consists  of  two  linear  FM  signals  mixed  with  white  Gaussian  noise.  In  terms  of  the 
variables  in  (2.27), 

^  r  "  n  2  n 

x(n)  =  Y^exp  } , 

i  •  1 

M  is  two,  Uj  =  8,  hj  =  94,  ^2  =  196,  ^2  =  30,  and  N  =  256.  With  a  record  length 
equal  to  one  second  and  a  sampling  rate  of  256  samples/second,  chirp  1  exhibits  a  slew 
rate  of  188  Hz/second,  starts  at  a  frequency  of  8/256,  and  finishes  at  196/256.  Likewise, 
chirp  2  has  a  slew  rate  of  60  Hz/second,  starts  at  a  frequency  of  196/256,  and  finishes  at  a 
frequency  equal  to  256/256.  The  parameter  u,  is  the  fractional  starting  frequency  of  the 
j*  chirp  times  N.  The  parameter  h,  is  equal  to  half  of  the  fractional  frequency  difference 
of  the  chiip  times  N,  or  simply  half  of  the  chirp’s  slew  rate. 

The  IHOMS  method  was  applied  first  by  fc .lowing  the  procedure  outlined  on  pages 
19  through  21.  Choosing  to  work  with  a  fourth  order  moment  slice  of  the  signal  fixes  the 
length  of  the  parameter  vector,  a,  to  four  (i.e.,  g  =  [  1,  Oj,  a2, 03] ).  GEN_A  is  an 
extrinsic  Matlab  function  in  Appendix  C  used  to  generate  a  desired  number  of  different 
parameter  vectors  having  elements  that  satisfy  (2.28)  and  (2.29).  Three  parameter  vectors 
are  formed  for  the  first  simulation.  The  remaining  steps  in  the  procedure  are  carried  out 
by  the  extrinsic  Matlab  function  ATH_IMS  in  Appendix  C.  One-hundred  lags  are  used  to 
compute  the  fourth-order  moment  slice  by  (2.30)  and  one  time-frequency  representation  is 
calculated  in  accordance  with  (2.31)  for  each  of  the  three  parameter  vectors.  The  value 
for  the  parameter  g  in  the  exponent  of  the  kernel  function  in  (2.31)  is  set  equal  to  0.0008 
because  the  same  value  was  found  to  give  acceptable  results  in  Example  1  of  [Ref.  6] .  The 
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final  time-frequency  representation  is  obtained  by  computing  the  magnitude  of  the  three 
parameter  vector  specific  surfaces.  Fifure  4.23(a)  shows  an  IHOMS  time-frequency 
representation  based  upon  a  fourth  order  moment  slice  and  three  parameter  vectors  for  a 
noise-free  signal  consisting  of  chirp  1  and  chirp  2.  The  two  chirps  are  indicated  by  the 
two  contour  lines  radiating  from  the  oi  igin.  The  remaining  lines  on  the  plot  are  cross¬ 
terms.  Figure  4.23(b)  is  a  plot  of  the  one-dimensional  function  0(0)  given  by 
(2.33).  This  function  represents  the  magnitude  of  the  time-frequency  surface  along  radial 
lines  passing  through  the  origin  of  the  surface.  The  first  maximum  in  Figure  4.23(b) 
occurs  at  0.09  radians  and  corresponds  to  the  b  parameter  of  chirp  2  through  (2.34); 

Applying  (2.34)  to  the  location  of  the  second  maximum,  0.28  radians,  an  estimate  of  the  b 
parameter  of  fre  first  chirp,  =  94.23,  is  obtained. 
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Figure  4.23:  Three  parameter  vector  IHOMS  representation  for  a  noise-free  signal. 
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Figxire  4.24  shows  the  three  parameter  vector  IHOMS  representation  for  the  same  two 
chirps  mixed  with  white  Gaussian  noise  at  an  input  SNR  of  -3dB.  The  chiips  cannot  be 
distinguished  from  the  background  noise. 


(a) 


(b) 


Figure  4.24;  Three  parameter  vector  IHOMS  representation  when  SNR  =  -3dB. 


Signal  detection  is  improved  by  summing  more  parameter  vector  specific  surfaces  in 
order  to  form  the  IHOMS  representation.  This  is  accomplished  by  creating  more 
parameter  vectors.  The  IHOMS  representation  in  Figure  4.25  uses  21  parameter  vectors 
on  the  two  chirp  signal  at  an  SNR  of  -3dB.  The  chirps  are  once  again  discernible.  Even 
with  21  parameter  vectors  the  detection  capability  degrades  significantly  as  the  SNR 
becomes  worse.  Whether  the  chirps  can  still  be  distinguished  in  Figure  4.26  when  the 
SNR  is  -5.5dB  is  questionable.  At  lower  SNRs  they  definitely  become  lost  in  the  noise. 

The  1  -  l/2Djps  -  spectrogram  estimates  of  the  double  chirp  signal  are  obtained  using 
ONE_HALF  and  SPECTRO  respectively.  A  rectangular  window  function,  a  32  point 
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'igure  4.25:  Twenty-one  parameter  vector  IHOMS  representation  when  SNR  =  -3dB. 
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Figure  4.26:  Twenty-one  parameter  vector  IHOMS  representation  when  SNR  =  -5.5dB 


data  window,  and  a  step  increment  of  one  are  the  specified  function  options  used  in  the 
simulations.  Figures  4.27,  4.28,  and  4.29  show  the  spectrogram  time-frequency 
representations  for  the  noise-free  signal,  an  SNR  of  -3dB,  and  an  SNR  of  -5.5dB, 
respectively.  Figures  4.30,  4.31,  and  4.32  display  the  1-1/2  Djpj  representations  for  the 
noise-free  case,  an  SNR  of  -3dB,  and  an  SNR  of  -5.5dB,  respectively.  Unlike  the 
IHOMS  representation,  both  1-1/2  Djp^  and  the  spectrogram  indicate  fractional  start  and 
finish  frequencies  for  each  chirp  in  addition  to  slew  rate.  For  a  given  noise  level,  the 
detection  capabilities  of  1-1/2  Djps,  the  spectrogram,  and  the  21  parameter  vector  IHOMS 
are  comparable.  However,  IHOMS  has  a  much  higher  computational  cost  than  either  of 
the  other  methods.  Generating  a  21  parameter  vector  IHOMS  surface  based  upon  a 
moment  slice  calculated  from  100  lags  requires  2,100  100-point  Fourier  transforms.  By 
comparison,  the  spectrogram  and  1-1/2  Djpj  require  only  256  32-point  Fourier  transforms 
to  generate  an  equivalent  time-frequency  representation. 


1200x 


(a) 


(b) 


Figure  4.27;  Spectrogram  of  noise-free  two  chirp  signal. 
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Figure  4.29:  Spectrogram  of  two  chirp  signal  with  SNR  =  -5.5dB 


Figure  4.32:  1-1/2  Djpg  representation  of  two  chirp  signal  with  SNR  =  -5.5dB. 


IHOMS  representations  indicate  only  the  slew  rate  of  chirps  and  not  their  start/stop 
frequencies.  As  a  consequence,  IHOMS  caimot  distinguish  between  two  or  more 
different  chirps  having  the  same  slew  rate.  Figure  4.33  displays  the  IHOMS 
representation  for  two  chirps.  The  first  chirp  starts  at  a  frequency  of  64/256  and  finishes 
at  192/256.  The  second  chirp  starts  at  128/256  and  fmishes  at  256/256.  They  both  have 
the  same  slew  rate  of  128Hz/second.  The  IHOMS  representation  detects  the  presence  of 
only  one  chirp  in  the  signal.  The  maximum  of  D(0)  occurs  at  0.19  radians  which 
corresponds  to  the  proper  b  parameter  of  64.  The  spectrogram  and  1-1/2  Djps 
representations  indicate  the  presence  of  both  chirps.  The  1-1/2  Djpj  time-frequency 
surface  of  this  signal  is  displayed  in  Figure  4.34. 
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The  IHOMS  algorithm  also  requires  a  more  intricate  sampling  scheme  than  does  the 
spectrogram  and  1-1/2  Djpg  methods.  The  received  signal  has  to  be  sampled  at  a  much 
higher  sampling  rate  in  order  to  obtain  the  data  samples  needed  to  form  the  moment  slices 
in  accordance  with  (2.30)  while  satisfying  (2.28)  and  (2.29). 
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V.  CONCLUSIONS 


A.  DISCUSSION  OF  RESULTS 

The  1*1/2  Dips  incthod  handles  both  signal  and  noise  power  differently  than  the 
spectrogram.  With  respect  to  signal  power,  the  maximum  processed  signal  power  is  the 
same  for  either  method.  But,  because  1-1/2  Djps  signal  power  is  modulated,  the  average 
signal  power  is  either  a  third  of  the  spectrogram’s  if  the  signal  is  real  valued,  or  half  the 
spectrogram’s  if  the  signal  is  complex  valued.  With  respect  to  noise  power,  the  methods 
exhibit  estimate  variances  that  depend  differently  on  the  size  of  the  data  window  and  the 
input  variance.  Whereas  the  spectrogram’s  variance  is  a  function  of  input  variance 
squared  and  window  length  squared,  1-1/2  Djp^’s  variance  depends  upon  input  variance 
cubed  and  window  length.  Consequently,  1- 1/2  Djps  is  mwe  resistant  to  input  noise  than 
the  spectrogram  if  the  noise  level  is  low  enough  and  the  window  length  is  large 
enough.  Conversely,  the  spectrogram  handles  noise  better  if  either  the  window  is  too 
small  or  the  input  noise  is  too  great.  When  the  effects  of  signal  and  r  Jse  are  combined, 
the  methods  are  compared  by  determining  how  many  noise  standard  deviations  the  average 
?’gnal  power  is  above  the  mean  noise  level.  This  measure  represents  the  method’s 
processing  gain.  1  - 1/2  Djps  has  a  larger  processing  gain  than  the  spectrogram  at  high  SNR 
for  any  reasonably  sized  data  window.  As  SNR  decreases,  the  processing  gain  advantage 
enjoyed  by  1-1/2  Dip,  diminishes.  At  a  sufficiently  low  SNR  value  that  depends  on  the 
size  of  the  data  window,  the  spectrogram  becomes  the  better  estimate.  As  window  length 
increases,  1-1/2  Djps’s  processing  gain  improves  relative  to  the  spectrogram’s. 

As  implemented  in  this  study,  the  spectrogram  and  1-1/2  are  superitM-  to  the 
bispectrum  at  detecting  multi-component  stationary  signals  in  white  Gaussian 
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noise.  However,  the  frequency  bin  averaging  scheme  used  in  the  simulations  is  not  as 
favorable  to  the  bispectrum  due  to  the  nature  of  its  frequency-frequency  representation. 
The  bispectrum  does  have  an  advantage  though  over  the  spectrogram  and  1-1/2  ^ips  in  its 
ability  to  distinguish  between  unrelated  signals,  harmonic  signals,  and  quadratically  phase 
coupled  signals.  The  only  limitation  is  that  quadratic  phase  coupling  is  not  detectable  if 
the  frequencies  of  the  coupled  signals  are  harmonics  of  the  same  fundamental  frequency. 

IHOMS  usefulness  is  limited  to  detection  of  linear  chirps  having  different  slew 
rates.  This  method  is  considerably  more  costly  than  the  spectrogram  and  1-1/2  Djp^  in 
terms  of  computations,  especially  for  noisy  signals. 

B,  SUGGESTIONS  FOR  FUTURE  STUDY 

Future  work  on  two  aspects  of  1-1/2  Djpg  might  prove  worthwhile.  The  first  is 
finding  an  optimum  window  and/or  smoothing  function.  Rectangular  windows, 
Hamming  windows,  and  boxcar  smoothing  have  been  used  in  ONE_HALF  to  date.  Since 
1-1/2  Djpj  is  a  degenerate  form  of  the  bispectrum,  a  window  function  satisfying  bispectral 
window  requirements  might  prove  to  be  better  window.  The  second  area  warranting 
further  study  involves  the  inherent  modulating  effect  of  the  1-1/2  Djps  algorithm.  The 
modulation  contains  useful  signal  information  if  it  could  be  extracted  from  the  1-1/2  Djp^ 
representation. 

The  logical  next  step  for  future  work  on  the  bispectruu  is  to  extend  it  to  a  time- 
frequency  representation.  Implementing  an  alternate  bispectral  calculation  method 
relying  on  polar  rasters  looks  like  a  promising  starting  point.  Evaluation  of  parametric 
methods  and  comparison  to  the  non-parametric  methods  of  computing  the  bispectrum 
might  prove  worthwhile.  In  addition,  better  averaging  schemes  for  the  bispectrum 
frequency-frequency  surface  could  be  developed. 
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APPENDIX  A:  VARIANCE  OF  THE  1-1/2  Djps  ESTIMATE 


The  1-1/2  Djps  spectnim  is  computed  using  (2.24).  The  variance  of  this  estimate  is 

derived  here  for  a  zero  mean  i.i.d.  Gaussian  input,  x  ~  N  (0,  a^) .  Results  are  obtained  for 

both  real  and  complex  valued  noise.  A  rectangular  window,  i.e.  w{k)=l  for  0  ^  1  ^  /V  -  1 , 
is  assumed  for  these  calculations. 

The  mean  of  the  1-1/2  Dips  estimate  is; 

£{1-1/2  t),ps}  =  £■{!-«(«)  V  (n-k)  +  X  (n)  {n  +  k) }  exp  (-jak) . 

^  km  0 

Since  the  expectation  operation  is  distributive. 


£{1-1/2  D.p,}  = 


N-l 

2  {E[\x(n)\^x*  (n-k)]  +Elx{n)\^x{n  +  k)]  }  expi-joik) . 


Let  term!  =  £{|x(n)|  V  (n-it) }  and  ter  m2  =  £  {  x  (n)  (n  +  A) }  .  Expand¬ 
ing  lerml  yields 

terml  =  £  {x^'  (/i)x(n)x’^  (n  -  A) }  . 


For  A  ^  0, 

terml  =  E  {x*  {n)xin)}  E  {x*  {n-k)} 
= 

X^X 

=  0, 


since  the  input  samples  are  independent  and  the  mean  of  the  input  sequence,  =  0. 
For  A  =  0, 

terml  =  £  {x"*  (n)x(n)x*  (n) }. 

Assuming  that  the  real  and  imaginary  parts  of  x(n)  are  independent. 
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Xfin),x^in)  ~NiO,  ^a^). 

then  term!  can  be  rewritten  as; 

terml=  E  {  [x^in)-jx-(n)]  [jc,(n)  +;>,  («)]  [jc,(n)  -;>,  («)]}  . 
which  simplifies  to, 

terml=  {£[j:J(n)]  +£  [jc,(/i)  ]£  [xj  (n)  ]}-;{£  [jrf  (/i)  ]  +£[x,(/i)]f  [-r?(«)  ] 
after  cross-multiplying  and  regrouping  the  real  and  imaginary  parts,  llie  terms  in  the  last 
expression  involving  the  product  of  a  first  moment  and  a  second  moment  are  zero  because 
the  first  moment  is  zero.  The  third  moment  quantities  in  the  same  expression  are  also 
zero  because  the  processes  are  distributed  normally  about  a  zero-mean.  The  final  result 
is  that  term!  is  zero  for  any  value  of  k.  The  same  is  true  for  ierml.  Consequently, 

£{  1-1/2  D/pj}  =  -  ^  [terml  +  term!}  exp  {-j(ji)k) 

=  0. 

The  second  moment  of  the  estimate  is: 

£{|l-l/2  =  £{1-1/2  bjpjin,w)  ■  1-1/2  (n,  to) } .  (A.l) 

where, 

—  (/n)  1  ,  *7  * 

1-1/2  Dips  («.  G))  =  2  ]C  in- m)  +  .x(n)rxin+m)}  exp  , 

m  •  0 

and 

1-1/2  Dip,  (n,co)  =  2  {  ix(n);^x(n-s) +!x(n)l^x  (n  +  s)}  expijcos)  . 

^sTo 

After  multiplying  the  last  two  equations  together,  collecting  terms,  and  distributing  the 
expectation  operator,  (A.l)  becomes; 
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, N-lN-l 
*•2  1 

£{|  1-1/2  D,pJ  ^  “  4  iC 

m  ■  Os  ■  0 

(A.2) 

where  terml  =  £[!jc(m)1^x*  (n  - /n)x(«  -  i)  ] , 

term!  =  £[iJt(/i)|\*  {n-m)x*  («  +  5)] 
term3  =  £  [jx(/i)|^jc  (w  +  m)  jf  («  -  5)  ] , 
and  termA  =  £  [|jt(/i);^jc(/j  +  ni)x*  (n  +  a)  )  . 

The  first  element  in  each  of  the  terms  above  represents  a  product  of  four  elements: 

ix(n)  ^  =  jc*  {n)x{n)x*  (n)x(n) . 

Terml,  through  termA  are  considered  separately  for  all  possible  m  and  s  index 
values .  T wo  different  mixed-moment  relationships  for  a  zero-mean  Gaussian  process  are 
needed.  For  a  real,  zero  mean  Gaussian  process,  the  following  relation  applies  [Ref.  7]; 

E{x^X2X^x^}  =  E{x^X2}E{x^x^)  ■¥E{x^x^)E{x2X^}  (A.3) 

+  E{x^x^)E{x2X^}  . 

For  a  complex  zero  mean  Gaussian  process,  a  slightly  different  relation  holds  [Ref.  11]: 
E{}t  XX2X* -iX^}  =  £  £{/ -^x^}  +£  {jf’^ijr4}  £  {/3X2}  .  (A.4) 

Terml  is  evaluated  first: 

1.  For  m  ^  0  and  5  0, 

terml  =  £  [j:”^  (n)  j:  («)jt(n)  ]  £  [x*  (n  -  m)x  («  -  a)  ]  .  (A.5) 

Applying  (A.3)  to  a  real  process,  terml  reduces  to: 

terml^  =  3  [£{jr(n)A:(n) }  J^£{x(n-m)x(n-s) }  (since  i.i.d.) 

=  3(ap  R^(s-m)  (since zero-mean) 
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=  3a^  ■  a^6  (s  -  m)  (since  ii.d.) 
=  3o^6  (s  -  m) . 


(A.6) 


For  a  complex  process,  (A.4)  is  employed  on  (A.5)  to  achieve  a  similar  result; 

terml^  =  2  [£  {jt(n)x*  (n) }  ]^E  {x*  {n-m)xin-s) } 

=  2ial)^R,{m-s) 

*  2o^  •  0^6  (w-s) 

-2a66(m-s).  (A.7) 

2.  If  m  =  0  and  s  0, 

terml^  =  £  {x*  (n)x  (n)x*  (n)x(/i)x*  (n)  }  £  {x{«  -  s) }  (A.8) 

=  0  (zero-mean  process  =❖£  {x(/j  -  j) }  =0) 

=  terml^. 

This  result  is  the  same  for  both  a  real  process  and  a  complex  process. 

3.  If  m  ^  0  and  s  =  0, 

terml^  =  £  {x*  (n)x(n)x*  (n)x(/i)x(n) }  £  {x*  («  -m) }  (A.9) 

=  0  (zero-mean  process  =>  £  {x’^  (n- m)  }  =  0) 
ss  terml^. 

4.  If  m  =  5  =  0, 

terml^  =  £  {x*  (n)x(n)x*  (n)x(n)x*  {n)x(n) } .  (A.IO) 

Assuming  that  the  real  and  imaginary  parts  of  the  process  are  independent 
implies  that, 

x^(/?),x,.(/i)  ~Af(0,  ^0^)  (A.11) 

Dropping  the  index  notation  in  (A.IO),  breaking  up  each  term  into  its  real  and 
imaginary  parts,  and  grouping  like  terms  produces; 

terml^  =  £{  (x^-yx,.)^(x,+yx,.)^}  , 
which,  after  multiplication  and  collection  of  terms,  becomes; 
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terml^  =  +  3£  {jcJ}  £  {jtJ}  +  3£  {jcJ}  £  {xf}  +£{xf}.  (A.12) 

The  second  order  moments  in  (A.12)  are  simply  the  variance  of  the  real  and 
imaginary  parts  of  the  process.  The  fourth  order  moments  are  computed  using 
(A.3).  The  sixth  order  moments  are  obtained  by  computing  the  sixth  deriva¬ 
tive  of  the  moment  generating  function  for  a  zero-mean  Normal  process  [Ref. 
17]: 


=  15o‘; 


'/-O 


-J: 

di^ 


£Xp 


ImO 


where  the  subscript  p  denotes  either  or  x,  since  their  variances  are  equal.  For 
a  complex  process,  in  terms  of  this  variance,  (A.12)  becomes: 


lerml'=  =  15o«  +  3(oJ)  (3a<)  +3(3o<)  (oj)  +  15oJ 
=  48ct«. 

From  (A.  11),  Using  this  relation,  (A.  13)  becomes: 

-  t 

termr  = 

=  6a!. 


If  the  process  is  real,  the  x,  terms  in  (A.12)  are  zero  and  x^  terms  become  x: 

terml^  =  £ {x^} . 

Evaluating  the  sixth  moment  using  the  moment  generating  function  yields: 

rerml^  =  15a®. 


(A.13) 


(A.14) 


(A.15) 


Term2  was  similarly  analyzed.  For  the  m ^ 0  and  s*Q  index  combination,  the 
second  moment  term  in  the  equivalent  (A.6)  and  (A.7)  is: 

£{x*  (/i-m)x(/i  +  5) }  =  a^6(-m-5) 

=  0, 
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because  the  delta  function  is  always  zero  for  non-negative  indices.  Like  ter  ml,  ter  m2  is 
zero  whenever  just  one  index  is  zero  because  the  process  is  zero  mean.  If  the  process  is 
real  and  both  indices  are  zero  then  term!  becomes  a  sixth  order  moment; 

terml^  =  E  {x (n) x  (n) x (n) x (n) x (n) x (n) } 

=  15aJ. 

If  the  process  is  complex,  four  terms  in  the  expectation  are  conjugated  while  two  are  not. 
The  result  is  a  final  expected  value  of  zero  since  the  process  is  zero  mean  Gaussian  with 
independent  real  and  imaginary  parts; 

terml^  =  E  [x*  (/i)x(n)x*  (n)x(n)x*  in)x*  (n) } 

=  £  {  -jx/) ^  (x^  +7jf,.)  2} 

=  0. 

Consequently,  term!  is  zero  for  all  index  combinations  if  the  process  is  complex  and  is 
non- zero  only  for  the  m  =  ^  =  0  condition  if  the  process  is  real. 

Evaluation  of  termS  produces  the  same  result  obtained  as  for  term!,  and  evaluation  of 
termA  produces  the  same  results  obtained  as  for  terml.  Table  A.l,  and  Table  A.2 
summarize  term  evaluation  results  for  the  real  and  complex  processes,  respectively. 


TABLE  A.1:  SUMMARY  OF  TERMS  FOR  REAL  PROCESSES 


index 

condition 

terml 

term2 

terms 

term4 

m*0  and 

3a^6(5-m) 

0 

0 

m  =  0  and 

0 

0 

0 

0 

m*0  and 

5  =  0 

0 

0 

0 

0 

m  =  s  -  0 

15a; 

15a; 

15aJ 

15a; 
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TABLE  A  SUMMARY  OF  TERMS  FOR  COMPLEX  PROCESSES 


index 

condition 

terml 

term2 

terms 

term4 

m*0  and 
s^O 

2a®6(5-m) 

0 

0 

m  =  0  and 
s*0 

0 

0 

0 

0 

m*0  and 

5  =  0 

0 

0 

0 

0 

0 

0 

6a® 

To  obtain  the  second  moment  of  the  1-1/2  Djpj  estimate  for  a  real  process,  the  terms 

shown  in  Table  A.  1  are  substituted  into  (A.2).  The  m  =  s  =  0  indices  are  taken  out  of 
the  summation  to  obtain: 


/V-l/V-l 

£{|  1-1/2  ^  {4(15aJ)  +  2(3ap5(s-m)exp(-J(o(m-s))} 

m  ■  Is**  1 

W-1 

=  ^  {  60a®  +  2  (3a®)  exp  (-jam) }  (sifting  theorem) 

m  ■  1 

=  ^{60aJ+(yv-l)6a5} 


54  +  6N.  2^3 
=  — ^  (op  ; 


(A.16) 


where  co  =  0  was  assumed  going  from  the  second  line  to  the  third  line  in  order  to  account 
for  the  maximum  variance  contribution  from  the  sum  term. 

Following  the  same  steps  for  a  complex  signal,  (A.2)  becomes: 

£  {|  1-1/2  }  COMPLEX  ~  [^■•■2]  (a^)  .  (A.17) 
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The  variance  of  the  1-1/2  Dips  estimate  is  given  by; 

var  {1-1/2  Dip,}  =  £{il-l/2  D^pJ^}  -  [£{1-1/2  D^p,}  ]^ 

Note  that  the  variance  is  simply  equal  to  the  second  moment  since  the  mean  of  the 
estimate  is  zero.  Therefore,  the  variance  of  the  1-1/2  Djp^  estimate  is  given  by  (A.  16)  for  a 
real  zero-mean  independent  Gaussian  process.  If  the  process  is  complex  zero-mean 
independent  Gaussian  with  independent  real  and  imaginary  parts,  the  estimate’s  variance 
is  given  by  (A.  17). 
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APPENDIX  B:  THEORETICAL  1-1/2  Djps  SIGNAL  POWER 


This  Appendix  first  shows  how  the  theoretical  1-1/2  Djpg  spectrum  is  calculated  for  an 
input  signal  consisting  of  a  single  real  valued  sinusoid.  The  computation  is  then  repeated 
for  a  complex  valued  sinusoid.  B  oth  calculations  are  made  using  a  rectangular  window  in 
(2.24). 

A.  REAL  VALUED  SINUSOID  CALCULATION 
The  real  valued  sinusoidal  input  signal  is 

x(n)  =  cos(0/i),  (B.l) 

where. 


0  =  2n^.  (B.2) 

is  the  digital  frequency  of  the  signal.  Substituting  (B.l)  into  (2.24)  using  a  rectangular 
window  yields: 

N-l 

1-1/2  Djpj  (n,  ®)  =  2  XI  [0  (”  -  ^)  ]  +  cos  [0  (n  +  /:)  ] }  exp  {-jwk) . 

kmO 


(B.3) 


After  applying  the  trigonometric  cosine  product  identity, 

cos  (A -B)  +  cos  (i4 +B)  =  2cosAcosB,  (B.4) 

and  pulling  the  cosine  cubed  term  out  of  the  summation  since  it  is  not  dependent  on  k, 
(B.3)  becomes: 


N-l 

1-1/2  Djpj  (n,  ©)  =  cos^  (0n)  ^  cos  (0/:)  exp  (-ycoA) .  (B.5) 

The  summation  portion  of  (B.5)  is  the  Fourier  transform  of  the  sinusoidal  signal. 
Evaluating  the  transform  at  ©  =  0  and  expressing  0  as  shown  in  (B.2),  the  Fourier 
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transform  can  be  expressed  in  terms  of  delta  functions; 

^cos  (In^k)  exp  i-j2n'^k)  =  ^  +5[a)- (jV-m)] } . 

*-0 

With  (B.2)  and  (B.6),  (B.5)  now  becomes: 

1-1/2  Djpj(«,©)  =  ^cos  (2n^«)  {6(oi)-m) +6[(o- (N-m)] }  . 


(B.6) 


(B.7) 


B.  COMPLEX  VALUED  SINUSOID  CALCULATION 
For  a  complex  valued  sinusoid  input. 


m 


(B.8) 


x(n)  =  exp(J2n-^n) 

=  exp  UQn) , 

the  1-1/2  Djps  estimate  calculated  by  (2.24)  using  a  rectangular  window  is; 

N-l 

1-1/2  Djpj.(n,co)  =*  2  E  \expijQn)\^  {exp[-jQin--k)]  +exp[jd(n  + k)]}  expi-jayk) 

(B.9) 

After  factoring  out  the  exponential  terms  that  are  not  a  function  of  k,  and  recognizing  that 
the  magnitude  squared  term  is  simply  equal  to  one,  (B.9)  reduces  to; 

1 

1-1/2  Djpj  (n,  co)  =  2  t  [-J0«] }  Yi  ^^P  ^^P  •  (®10) 

i-O 


With  the  Euler  identity. 


1 


2  i^^P  L/0”]  +  ^^P  }  =  cos  (dn) , 


(B.ll) 


and  (B.2),  (B.IO)  becomes; 

1-1/2  Djpj.(/J,©)  =  cos(27i^n)  exp(j2n’j^k)exp{-j(iik) 


m 


N-l 


m , 


kmO 


(B.12) 


The  summation  term  in  (B.12)  is  the  Fourier  transform  of  a  complex  sinusoid  which  is 
expressible  in  terms  of  a  delta  function  yielding; 


1-1/2  (d)  = /Vcos(2n^n)6(co-m) . 

ips  yV 


(B.13) 
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APPENDIX  C:  COMPUTER  PROGRAMS 


The  four  extrinsic  Matlab  functions  contained  in  this  Appendix  are  written  for  Matlab 
version  4.0.  If  the  functions  are  used  with  earlier  versions  of  Matlab,  contour  plots  and 
mesh  plots  will  not  be  oriented  correctly. 
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%function  [P,freqindex,timeindex]  =  one_half(data,wintype,winlen,step); 

%This  function  will  calculate  the  1  1/2  Dips  spectral  surface.  The  1  1/2  Dips 
%surface  characteristics  are  determined  by  the  selection  of  window 
%type  (wintype),  window  length  (winlen)  and  the  distance  that  the  window 
%is  moved  through  the  data  sequence  (step).  The  magnitude  of  the  positive 
%half  of  the  1  1/2  D  spectral  plane  is  returned  in  the  “P”  matrix. 

%Ouq)uts  “timeindex”  and  “freqindex”  aid  in  axis  labeling. 

%The  inputs  are: 

%data  -  The  input  observations  vector.  The  length  should  be  a  power  of  2. 
%wintype;  ‘0’  Rectangular  Window 
%  ‘1’  Hamming  Window 

%winlen  -  The  desired  width  of  the  window,  normally  half  of  the  input 
%  length. 

%step  -  Time  step  desired,  can  be  ‘1’  or  a  multiple  of  ‘2’ 

%prepared  by  Karen  A.  Hagerman,  06  May  1992. 

%modified  by  Jeff  McAloon,  01  June  1993. 


function  [P,freqindex,timeindex,svect]=one_half(data,wintype, winlen, step) 

[datarows,datacolumns]=size(data); 
if  datarows-^l 
data=data.’; 


end 

siglen=length(data); 
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if  wintype=0 

win=ones(winlen- 1,1); 
elseif  wintypc==l 
win=hamming(winlen- 1 ); 
end 

W=[wm(winlen/2:- 1:1)]; 
x=[zeros(l,winlen)  data  zeros(l,winlen)].’; 
p=reros(siglen/step,winlen); 


index=l; 

for  n=winlen+l  :step:siglen+winlen-step+l 
xt=x; 

U=max((winien+ 1  ),(n-(winlsn/2- 1 ))); 

ul=min((n+(winlen/2- 1  )),(siglen+winlen)); 

ctp=xt(ll:ul); 

mtp=mean(ctp); 

xt(ll:ul)=ctp-mtp; 

xdt=xt((n-(winlen/2- 1 )) :  (n+(winlen/2- 1 ))); 
Xn=[abs(xt(n))''2;abs(xt(n))''2]; 
Xm=[conj(xt(n:-l;n-(winlen/2-l)))  xt(n:n+(winlen/2-l))] 
product=((Xm*Xn)*W).’; 
product=[product  0  conj(product(winlen/2:-l:2))]; 
pt=0.5*fftshift(fft(product)); 
p(index,:)=pt; 
index=index+l; 
end 

P=abs(p( :  ,winlen/2+ 1 :  winlen)); 
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[prow,pcoluinn]=size(P); 
freqmdex=[0:pcolumn- 1]; 
timcindex=[  1  ;prow] ; 


%  function  [P.freqaxis.timeaxis]  =  spectro(data,wintype,winlen); 

% 

%  This  function  calculates  spectrogram  of  the  supplied  sequence, 

%  “data”.  The  user  must  specify; 

%  “wintype”  -  “0”  for  a  rectangular  window. 

%  “  1  ”  for  a  Hamming  window. 

%  “winlen”  -  The  desired  length  of  the  data  window. 

%  The  time  step  is  fixed  at  one  and  the  spectrogram  is  calculated  with 
%  non-normalized  periodograms.  The  time-frequency  surface  is  returned 
%  in  the  “P”  matrix.  The  columns  of  “P”  are  the  frequency  bins  while  the 
%  rows  are  the  time  steps.  Time-frequency  axis  labeling  vectors, 

%  “freqaxis”  and  “timeaxis”  are  also  returned  to  aid  in  the  plotting 
%  of  results. 

%  prepared  by  Jeff  McAloon,  01  June  1993. 


function  [pow,freqaxis,timeaxis]  =  spectro(data,wintype,winlen) 

x=data; 

xlen=length(x); 

N=winlen; 

[row,  col]=size(x); 
if  row  >  1 


end 

x=[zeros(l,N/2)  x  zeros(l,N/2)]; 
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if  wintype  ==  1 

win=hamming(N)  ’ ; 
elseif  wintype  =  0 
win=ones(l,N); 
end 

for  ind=l:xlen 

xseg=x(ind;(N-  1+ind)); 
xseg=xseg-mean(xseg) ; 
SP=(abs(fft(xseg.*win,N))).'^2 
POW(ind.0=fftshift(SP); 
end 

pcw=P0W(:J^/2+l;N); 
tinieaxis=0:xlen-l; 
freqaxis=0:  N/2- 1 ; 


%  function  [bis,  baxis]  =  indbis(data,datseg,lag,tlen,win) 

%  This  function  computes  the  bispectrum  of  real  and  complex  valued  data 
%  sequences  by  the  indirect  method.  The  tri-correlation  sequence  is 
%  computed  for  a  specified  number  of  segments.  These  segment  sequences 
%  are  then  averaged  to  obtain  the  final  tri-correlation  sequence.  The 
%  three  available  tri-correlation  window  options  are  the  unit  hexagonal 
%  window,  Parzen’s  window,  and  the  optimum  window  (also  known  as  Sazaki’s 
%  minimum  bias  window).  The  required  function  call  arguments  are: 

%  “data”  ->  input  data  sequence  vector. 

%  “datseg”  ->  number  of  samples  to  be  used  in  each  segment. 

%  “lag”  ->  number  of  lags  to  be  used  in  the  computation  of  the 
%  tri-correlation  sequence. 

%  “tlen”  ->  square  dimension  of  the  two  dimensional  FFT  to  be  used 
%  on  the  final  form  of  the  tri-correlation  sequence. 

%  “win”  ->  “0”  for  unit  hexagonal  window 
%  “1”  for  optimum  window 

%  “2”  for  Parzen’s  window. 

%  The  function  returns  the  complex  bispectrum  in  the  “tlen-by-tlen”  sized 
%  matrix  “bis”.  A  vector  “baxis”  is  returned  that  can  be  used  to  label 
%  both  bispectrum  frequency  axes. 

%  prepared  by  Jeff  McAIoon,  01  Jun  1993. 

function  [bis,  baxis]  =  indbis(data,datseg,lag,tlen,win) 

%  Initialize  parameters  and  reorient  data  if  necessary: 
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N=length(data); 

Ns=fix(N/datseg); 

[rd,  cd]=size(data); 
if  cd  >  1 
x=data.’; 
else 
x=data; 
end 

x=[0;x-mean(x)]; 

Rt=zeros(2*lag+l); 
rx=zeros(2*lag+ 1 ) ; 

%  Compute  tri-correlation  sequence; 
for  i=l:Ns 
for  jslidatseg 

y(j)=x(j+(i- 1  )*datseg); 
end 

y=y-mean(y); 
for  m=l:(2*lag+l) 
for  n=l:(2*lag+l) 

s  1  =max([  1 (m-  (lag+ 1  )),-(n-(lag+ 1 ))] ); 
s2=min([lag,lag-(m-(lag+l)),iag-(n-(lag+l))]); 

p=0; 

fork=sl;s2 

r=r+conj(y(k+l))*y(k+l+(m-(lag-i-l)))*y(k+l+(n-(lag+l))); 

end 

rx(n4n)=r; 

end 
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end 

Rt=Rt+rx; 

«=□; 

end 

R=flipud(Rt/Ns); 

R=Rt/Ns; 

%  Determine  tri-correlation  window  function: 

if  win  =  0  %  unit  hexagonal  window 

W=ones(si2e(R)); 

elseif  win  =1  %  optimum  window 

itm=0; 

for  m=-lag:lag 
itm=itm+l; 

wm=abs(sin(pi*m/lag))/pi  +  (l-abs(m/lag))*cos(pi*m/lag); 
itn=0; 

for  ns-lag:lag 
itn=itn+l; 

wn=abs(sin(pi*n/lag))/pi  +  (l-abs(n/lag))*cos(pi*n/lag); 
wmn=abs(sin(pi*(m-n)/lag))/pi  +  (l-abs((m-n)/lag))*cos(pi*(m-n)/lag); 
W(''tm,itn)=wm*wn*wmn; 
end 
end 

elseif  win  =  2  %  Parzen’s  window 

itm=0; 

for  m=-lag:lag 
itm=itm+l; 
if  abs(m)  <=  lag/2 
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wm=l-6*(abs(m)/lag)''2+6*(abs(m)/lag)''3; 

else 

wm=2*(  1  -abs(m)/lag)''3; 
end 
itn=0; 

for  n=-lag:lag 
itn-itn+1; 
if  abs(n)  <=  lag/2 

wn=l-6*(abs(n)/lag)''2+6*(abs(n)/lag)''3; 
wmn=  1  -6*(abs(m-n)/lag)''2+6*(abs(m-n)/lag)''3; 
else 

wn=2*(  1  -abs(n)/Iag)''3; 
winn=2*(l-abs(m-n)/lag)'^3; 
end 

W(itm,itn)=wm*  wn*  wmn; 
end 
end 
end 

%  Compute  bispectrum. 
bis=fftshift(fft2(R.*W,tlen,tlen)); 
baxis=(-tlen/2);  (tlen/2- 1 ); 
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%  function  Amtrx  =  gen_a(al_eleni); 

%  This  function  computes  the  parameter  vectors  needed  to  compute  a  fourth 
%  order  moment  slice  in  the  IHOMS  time-frequency  method  implemented  by 
%  the  extrinsic  matlab  function  ATH_IMS.M.  The  first  element,  al,  of 
%  each  IHOMS  parameter  vector  is  specified  in  the  function  call  as  the 
%  vector,  “al_elem”.  The  length  of  “al_elem”  is  the  number  of  distinct 
%  parameter  vector  sets  (i.e.one  set  would  be  [1  al  a2  a3])  returned  as 
%  columns  in  the  matrix  “Amtrx”.  This  also  equals  the  number  of  time- 
%  frequency  surfaces  summed  by  IHOMS  to  form  the  fmal  time-frequency 
%  representation. 

%  prepared  by  Jeff  McAloon  01  June  1993. 

function  Amtrx  =  gen_a(al_elem) 

%  Verify  input  vector  a  column  vector; 

[rin,  cin]=size(al_elem); 
ifcin>  1 

alv=al_elem.’; 

end 

%  Iteratively  solve  for  parameter  vectors: 

A=[]; 

for  p=l:length(alv) 
a3v=-5;l:5; 

a2a=100;  a2b=l;  i=l;  md=99; 
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while  (abs(a2b-a2a)  >=  0.1)  &  (i  <=  length(a3v)) 
a2a=l-alv(p)-a3v(i); 
a2b=Teal(sqrt(l  -al  v(py^2-a3v(i)'^2)); 
if  abs(a2a-a2b)  <=  md 
in(lmin=i; 
end 

md=min(md,abs(a2a-a2b)); 

id=?nax(2,indmin); 

i=i+l; 

end 

a3v=a3v(id- 1):0.0001  ;a3v(id); 
a2a=l;  a2bs2;  i=0; 

while  (abs(a2b-a2a)  >=  0.0001)  &  (i  <=  (length(a3v)-l)) 
i=i+l; 

a2a=l-alv(p)-a3v(i); 
a2b=sqrt(l-al  v(p)A2-a3v(i)^2); 
end 

al*alv(p);  a2=a2a;  a3=a3v(i); 
a=[l  al  a2  a3]; 

A=[A,a.’]; 

end 

Amtrx=A; 
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%  function  surf  =  ath_ims(M,u,b,N,SNR.A,g) 


%  This  function  computes  the  IHOMS  time-frequency  representation  of 
%  multi-component  chirp  signals  in  white  Gaussian  noise.  Input  arguments 
%  are: 

%  M  ->  number  of  chirp  components. 

%  u  ->  vector  of  chirp  starting  frequencies. 

%  b  ->  vector  of  halved  chirp  slew  rates. 

%  N  ->  data  record  length. 

%  SNR  ->  desired  SNR.  If  noise  free  signal  desired  enter  “99”. 

%  A  ->  matrix  whose  columns  are  the  “a”  parameter  vectors  needed  to 
%  form  the  moment  slice.  See  extrinsic  Matlab  function  GEN_A 
%  generate  this  matrix. 

%  g  ->  kernel  fimction  parameter. 

%  The  moment  slice  is  computed  using  100  lags.  The  IHOMS  surface  is 
%  plotted  along  with  a  1-D  radial  maxima  plot  that  aids  in  locating  and 
%  characterizing  the  chirps  present  in  the  signal.  In  addition,  the 
%  function  returns  a  100  X  100  matrix,  “surf’.  The  columns  are  time/lag 
%  and  the  rows  are  frequency  bins. 

%  prepared  by  Jeff  McAloon,  01  June  1993. 

function  surf  =  ath_ims(M,u,b,N,SNR,A,g) 

%  Initialize  parameters: 
randn(‘seed’,0) 

if  SNR  =99  %  noise  free 
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GG=0; 

else 

GG=1; 

end 

Z==sqrt(10^(GG*SNR/10)); 

maxlag=100; 

j=sqrt(-l); 

[rA,cA]=sizc(A); 

SURF=zeros(  100); 

%  IHOMS  algorithm: 
for  k=l:cA 
a=A(;.k); 
kn=k 
1^=0; 

for  n=0:maxlag*l 
for  lag=0:maxlag-l 

no=GG*ran(in(  1  ,length(a)); 

term(l)=0; 

for  p=l:M 

term(l)=term(l)+Z*exp(j*2*pi*(((lag)/N)*u(p)+(((lag)/N)''2)*b(p))); 

end 

term(  1  )=term(  1  )+no(  1 ) ; 
for  i=2:length(a) 
args=n+a(i)*lag; 
if  arg<0 
term(i)=0; 
else 
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tcnn(i)=0; 

forq=l:M 

term(i)=term(i)+Z*exp(j*2*pi*(((arg)/N)*u(q)+(((arg)/N)'^2)*b(q))); 

end 

term(i)=tenn(i)+no(i); 

end 

end 

kem=exp(-g*lag'^2); 

m(lag+l)=kem*conj(temi(l))*term(2)*term(3)*tenn(4); 

end 

TF(n+l.:)=fft(m); 

end 

SURF=SURF+TF; 

end 

%  Orient  matrix  for  MATLAB  4.0  contour  plot: 
surf=rot90(flipud(abs(SURF)),- 1 ); 

%  Plot  IHOMS  surface; 

t=0:maxlag-l; 

f=0;maxlag-l; 

subplot(121),contour(t,f,surf,2),title(‘IHOMS  surface’), ... 
xlabel(‘time/lag’),ylabel(‘frequency’) 

%  IHOMS  1_D  search  routine  and  plot: 

theta=0:0.01:pi/2; 

for  in=l:length(theta) 

D(in)=0; 


109 


for  r=30:90 


t=  1  +r*cos(theta(in)); 
f=  1  +r*sin(theta(in)); 

D(in)=sD(in)+surf(f,t); 

end 

end 

subplot(122),plot(theta,D)^abel(‘theta  in  radians’). 
ylabel(‘D(theta)’),title(‘l-D  radial  maxima  plot’) 
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